ANEOS SILICA MODEL DOCUMENTATION NOTEBOOK¶

ANEOS FUSED SILICA v0.2¶

This Jupyter notebook provides documentation for a preliminary ANEOS model for fused SiO$_2$.

Reference: Sarah T. Stewart and Kaitlyn M. Amodeo (2024). ANEOS Equation of State Model for Fused Silica. Zenodo. http://doi.org/10.5281/zenodo.XXX (https://github.com/ststewart/aneos-fused-silica/).

When this model EOS is used, cite this doi-tagged Zenodo GitHub repository. Future updates to this EOS will be archived in Zenodo doi-tagged version-controlled releases.

Sarah T. Stewart sstewa56@asu.edu October 17, 2024

VERSION INFORMATION
EOS Version: Fused_Silica-ANEOSv0.2G1
SLVT# translates to solid-liquid-vapor-thermal model adjustment-version number
G1 = the gridded EOS tables used the gridflag=1 definition for the grid points
SLVTv0.2G1 Release Date: 20241017
GitHub: https://github.com/ststewart/aneos-fused-silica/

This is a solid-liquid-vapor ANEOS model with ionization for SiO$_2$ glass. This model EOS uses the heat capacity adjustment factor described in Stewart et al. (AIP Conference Proceedings, 2020) and Stewart, ANEOS Code Modification: Thermal model adjustment parameter, https://github.com/ststewart/aneos-forsterite-2019/EOS-docs/, 2019.

This Jupyter notebook was used to develop and document ANEOS input parameters that are optimized to try to span very low pressures and temperatures to very high pressures and temperatures.

USAGE
This notebook is designed to be used in two modes:

  • Development mode: ANEOS is called in the local directory to generate new model output. The user can varying the model input parameters and tabular gridding schemes. This notebook produces plots of the equation of state surface and comparisons to experimental data and ab initio calculations.
  • Release mode: the notebook interacts with previously created ANEOS model output. This notebook provides detailed documentation of a version-controlled release of the model EOS and enables users to interact with the EOS surface and compare the model to data.

This notebook is distributed in release mode. Users who wish to use the development mode can obtain the ANEOS source code from: Thompson, S. L., Lauson, H. S., Melosh, H. J., Collins, G. S., & Stewart, S. T. (2019, November 1). M-ANEOS (Version 1.0). Zenodo. http://doi.org/10.5281/zenodo.3525030.

Development mode:
The notebook calls a local aneos executable which expects a formatted ANEOS.INPUT ascii file and an ascii gridtable.txt definition file that is generated by this notebook. ANEOS produces an ANEOS.OUTPUT ascii file and gridded ascii SESAME-format EOS tables.
The user edits (1) the input parameters in ANEOS.INPUT (separately from this notebook) and (2) the output EOS table parameters in the notebook cell below.
This notebook reads some parameters from ANEOS.INPUT and assumes that the variables are spaced by 10 columns.
When the notebook is run, the code issues a system call to run the aneos program in the local directory and reads in results from the ANEOS.OUTPUT, NEW-SESAME-STD.TXT, and NEW-SESAME-EXT.TXT files to make the plots below and to generate the GADGET format EOS table and other miscellaneous files.

Release mode:
This notebook reads in the provided ANEOS.INPUT, EOS tables, and other distribution data to generate the plots below.

OUTPUT/PROVIDED DATA FILES
ANEOS.OUTPUT: Plain text file with an overview of the calculated equation of state. Look at this file for any errors in the EOS calculation.

NEW-SESAME-STD.TXT: Standard length Sandia-style SESAME file with 201 table and 301 table (density, temperature, pressure, sp. internal energy, Helmholtz free energy). 301 table units: g/cm$^3$, K, GPa, MJ/kg, MJ/kg.

NEW-SESAME-EXT.TXT: SESAME-style table with extra variables from ANEOS. Contains the standard 201 table and non-standard 301-extra-variables EOS table. The 301 table has: density grid values, temperature grid values, sp. entropy(T,rho), sound speed(T,rho), sp. heat capacity(T,rho), KPA flag(T,rho). 2-D arrays list all densities, looping over each temperature. 301 table units: g/cm$^3$, K, MJ/K/kg, cm/s, MJ/K/kg, integer flag, integer flag. The KPA flag is an ANEOS output with phase information.

NEW-SESAME-HUG.TXT: Ascii table with the reference state Hugoniot.

NEW-SESAME-MDQ.TXT: Ascii 301-style SESAME table with the MDQ rating.

NEW-SESAME-IEP-BILINEAR.TXT: Standard length Sandia-style SESAME file with 201 table and 301-style table containing the IEP values for bilinear interpolation. Format: 1-D interpolated density points, 1-D interpolated temperature points, 2-D IEP arrays for pressure, sp. internal energy, Helmholtz free energy, sound speed). IEP units are ratios to the original calculated ANEOS point values. Located in the interpolation subdirectory.

NEW-SESAME-STD-NOTENSION.TXT: Standard length Sandia-style SESAME file with 201 table and 301 table (density, temperature, pressure, sp. internal energy, Helmholtz free energy) where the ANEOS tension region is removed and replaced with the solid-vapor coexistence region. 301 table units: g/cm$^3$, K, GPa, MJ/kg, MJ/kg.

NEW-SESAME-EXT-NOTENSION.TXT: SESAME-style table with extra variables from ANEOS. Contains the standard 201 table and non-standard 301-extra-variables EOS table where the ANEOS tension region is removed and replaced with the solid-vapor coexistence region. The 301 table has: density grid values, temperature grid values, sp. entropy(T,rho), sound speed(T,rho), sp. heat capacity(T,rho), KPA flag(T,rho). 2-D arrays list all densities, looping over each temperature. 301 table units: g/cm$^3$, K, MJ/K/kg, cm/s, MJ/K/kg, integer flag, integer flag. The KPA flag is an ANEOS output with phase information.

NEW-SESAME-GADGETINIT-NOTENSION.TXT: Standard length Sandia-style SESAME file with 201 table and 301 table (density, temperature, pressure, sp. internal energy, specific entropy) where the ANEOS tension region is removed and replaced with the solid-vapor coexistence region. 301 table units: g/cm$^3$, K, GPa, MJ/kg, MJ/K/kg.

NEW-GADGET2-STD-NOTENSION.TXT: Standard GADGET2 EOS table generated by interpolating the ANEOS table with no tension. Format: number of density points, number of sp. entropy points, density grid values, sp. entropy grid values, 2-D arrays of pressure, temperature, specific internal energy, sound speed. 2-D array values list all densities, looping over each value for specific entropy. Units: g/cm$^3$, erg/K/g, dynes/cm$^2$, K, erg/g, cm/s.

NEW-GADGET2-EXT-NOTENSION.TXT: Extended variables GADGET2 EOS table generated by interpolating the ANEOS table with no tension. Format: number of density points, number of sp. entropy points, density grid values, sp. entropy grid values, 2-D arrays of Helmholtz free energy, KPA flag, MDQ flag. 2-D array values list all densities, looping over each value for specific entropy. Units: g/cm$^3$, erg/K/g, erg/g, integer flag, integer flag.

The SESAME tabular equation of state formats are described in Lyon and Johnson (1992).

APPLICATIONS AND LIMITATIONS

The reference state is the solid fused silica at STP. In ANEOS, the liquid model is a modification of the solid model. As a compromise to span the whole phase diagram, the solid does not represent the true polymorphs of SiO$_2$. The solid phases are simplified into two phases: a low-pressure phase and a denser, high-pressure phase. The heat capacity adjustment parameter is adjusted to match the estimated heat capacity of the liquid as determined from shock temperature data. The calculated shock temperatures provide a good fit to the data. The model Debye tempereature is not the true value for the solid; the Debye temperature was fitted to attain the desired specific entropy at boiling at 1 bar. The temperatures in the liquid field are much improved with the adjusted heat capacity compared to previous ANEOS models using the Dulong-Petit limit.

This high-pressure phase transition introduces artificial discontinuities in the liquid field and along the melt curve. This model becomes increasingly erroneous for shock pressures above about 800-1000 GPa as the curvature on the Hugoniot and the Grueneisen parameter are not correct at extremely high densities. The functional form for the Gruneisen parameter in ANEOS has difficulty matching the properties of silicate liquids, and this model does not approach the Thomas-Fermi limit.

The experimental constraints on the melt curve for SiO$_2$ are poor. This preliminary model is tuned to agree at low pressures and is likely unreliable at high pressures.

The calculated vapor curve is skewed with a high temperature and low density, although it is close to the 1 bar values.

ANEOS NOTES
This work uses v1.0 of M-ANEOS released on Zenodo: Thompson, S. L., Lauson, H. S., Melosh, H. J., Collins, G. S., & Stewart, S. T. (2019). M-ANEOS (1.0). Zenodo. https://doi.org/10.5281/zenodo.3525030 This version included Melosh's (2007) treatment for molecular gas and the capability to include a melt curve and solid-solid/liquid-liquid transition (Collins & Melosh LPSC 2014). This version includes an adjustment to the Debye model for the thermal term in the Helmholtz free energy to approach a user-defined heat capacity at high temperatures. The multiplicative factor $f_{cv}$ is entered in input value V44, and the high-temperature heat capacity is $3f_{cv}Nk$.

The ANEOSTEST.f routine was modified to output tabulated EOS. Note that the current version of this function sets positive pressures smaller than 1.E-30 GPa equal to 1.E-30 GPa.
ANEOS2.f was modified to increase the number of points tabulated on the melt curve in the ANEOS.OUTPUT file and to gather the variables for the heat capacity modification.
ANHUG.f was modified to output more Hugoniot points.
ANEOS1.f and ANEOS2.f were modified to increase the high temperature limit for the heat capacity (Stewart et al., AIP Conference Proceedings, 2020).

CORRECTIONS AND IMPROVEMENTS
Please send corrections to STS and any requests for data to include in the model-data comparison plots.

In [1]:
# defines class objects for extEOStable and functions to read in ANEOS input and output
# READ THIS CODE IF YOU WANT TO UNDERSTAND THE OBJECT STRUCTURES

%run import-modules
Platform:  Darwin Kernel Version 24.0.0: Tue Sep 24 23:37:25 PDT 2024; root:xnu-11215.1.12~1/RELEASE_ARM64_T6030
python version:  3.11.9
matplotlib version:  3.9.1
hvplot version:  0.10.0
numpy version:  2.0.1
pandas version:  2.2.2

In [2]:
# THIS CELL DEFINES FUNCTIONS AND CUSTOMIZATIONS
#
plt.rcParams["figure.figsize"] = (20,10)
plt.rcParams["patch.force_edgecolor"] = True
#
# define index function (must be a more efficient way to do this...)
get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x == y]
#
# this prevents some error messages for inline plotting in a jupyter notebook
%matplotlib inline 
#
def make_error_boxes(ax, xdata, ydata, xerror, yerror, facecolor='r',
                     edgecolor='None', alpha=0.5):

    # Loop over data points; create box from errors at each point
    errorboxes = [Rectangle((x - xe[0], y - ye[0]), xe.sum(), ye.sum())
                  for x, y, xe, ye in zip(xdata, ydata, xerror.T, yerror.T)]

    # Create patch collection with specified colour/alpha
    pc = PatchCollection(errorboxes, facecolor=facecolor, alpha=alpha,
                         edgecolor=edgecolor)

    # Add collection to axes
    ax.add_collection(pc)

    # Plot errorbars
    artists = ax.errorbar(xdata, ydata, xerr=xerror, yerr=yerror,
                          fmt='None', ecolor='k')

    return artists
In [3]:
# Remember the versions used to run python code
# copy this cell into your own notebooks
# add additional modules to the list in the f.write statements

# Record the verions information for these calculation for posterity
from platform import python_version
from matplotlib import __version__ as mplversion
from numpy import __version__ as npversion
from scipy import __version__ as scipyversion
from pandas import __version__ as pdversion
from eos_table import __version__ as etabversion

# Get date and time values in ISO 8601 format
# importing datetime and time module
# Calling the isoformat() function over the current time including time zone:
import datetime
import zoneinfo
DateTime_in_ISOFormat = datetime.datetime.now(zoneinfo.ZoneInfo('America/Los_Angeles')).isoformat('T', 'seconds') #

versionfilename = './python-versions.txt'
with open(versionfilename, 'a') as f:
    f.write('---------------------------\n')
    f.write(DateTime_in_ISOFormat+'\n')
    f.write('python '+python_version()+'\n')
    f.write('numpy '+npversion+'\n')
    f.write('scipy '+scipyversion+'\n')
    f.write('matplotlib '+mplversion+'\n')
    f.write('pandas '+pdversion+'\n')
    f.write('eos_table '+etabversion+'\n')
    f.write('---------------------------\n')

# Printing Today's date and time in ISO format
print(f"Version history file {versionfilename}")
print(f"Date and time {DateTime_in_ISOFormat}")
print('python '+python_version())
print('numpy '+npversion)
print('scipy '+scipyversion)
print('matplotlib '+mplversion)
print('pandas '+pdversion)
print('eos_table '+etabversion)
Version history file ./python-versions.txt
Date and time 2024-10-18T08:51:18-07:00
python 3.11.9
numpy 2.0.1
scipy 1.14.0
matplotlib 3.9.1
pandas 2.2.2
eos_table 1.1.4

USER INPUTS FOR SESAME AND GADGET2 TABLE CONSTRUCTION¶

If the code cell below is hidden, use the button above to reveal the cell.

In development mode, the user must input:

  1. Header information for the SESAME table.
  2. Temperature, density and entropy grid points.

The following code cell also includes the development mode flag and option to skip construction of a GADGET2 table if it is not needed.

In [4]:
# -------------------------------------------------------------
# USER INPUTS
#-------------------------------------------------------------
# NOTEBOOK MODE
#-------------------------------------------------------------
#
DEVMODE = 0 # 1=development mode; 0=release mode
#
#-------------------------------------------------------------
# OPTIONS TO REMOVE TENSION REGION AND CONSTRUCT A GADGET2 TABLE
#-------------------------------------------------------------
#
REMOVETENSIONFLAG=1  # 0- keep tension region; 1- remove tension region
MAKEGADGETTABLE=1  # 1=make table; 0=do not make gadget2 table
#
#-------------------------------------------------------------
# USER MUST SET THESE VALUES TO DEFINE THE TABULATED EOS GRIDS
#-------------------------------------------------------------
#
# ====>>>>>> YOU NEED TO MAKE SURE THESE VALUES MATCH ANEOS.INPUT  <<<<=====
MATERIALNAME = 'silica'
MODELNAME = 'Fused-silica-ANEOS-SLVTv0.2G1'
# Header information must all be compatible with float format
MATID = 1.0        # MATID number
DATE = 220227.     # Date as a single 6-digit number YYMMDD
VERSION = 0.1      # ANEOS Parameters Version number
FMN = 30.          # Formula weight in atomic numbers for silica
FMW = 60.083       # Formula molecular weight (g/mol)
# The following define the default initial state for material in the 201 table
R0REF   = 2.20     # g/cm3 *** R0REF is inserted into the density array
K0REF   = 4.0E11   # dynes/cm2
T0REF   = 298.     # K -- *** T0REF is inserted into the temperature array
P0REF   = 1.E6     # dynes/cm2 -- this defines the principal Hugoniot calculated below
#
# generate output files for other ANEOS models
modelflag = 0   # 0 - don't change anything above; 1 - Melosh MAPS 2007
#
if modelflag == 1:
    MODELNAME = 'Silica_Melosh_2007'
    # Header information must all be compatible with float format
    MATID = 1.0        # MATID number
    DATE = 070513.     # Date as a single 6-digit number YYMMDD
    VERSION = 1.0      # ANEOS Parameters Version number
    FMN = 30.          # Formula weight in atomic numbers for silica
    FMW = 60.083       # Formula molecular weight (g/mol)
    # The following define the default initial state for material in the 201 table
    R0REF   = 2.65     # g/cm3 *** R0REF is inserted into the density array
    K0REF   = 3.58E11  # dynes/cm2
    T0REF   = 298.     # K -- *** T0REF is inserted into the temperature array
    P0REF   = 1.E6     # dynes/cm2 -- this defines the principal Hugoniot calculated below
#
#-------------------------------------------------------------
# EXAMPLE EOS GRID ARRAYS
#-------------------------------------------------------------
#
# Life is easier if the grid arrays are even in length (after adding T0REF and R0REF below)
#
# Flag for easily switching between different table grids during model development
# 0: low res for testing
# 1: high res for production
gridflag = 11
#
#------------------------------------------------------
if gridflag == 99: # low res for testing
    gridtvals  = [0, 10000., 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [30., 30.,20.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    gridrvals  = [1.e-20, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [20, 10, 6, 30, 20] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    gridsvals  = [1.e-4,.25] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [0] # 1 for linear, 0 for log in each segment
    gridsnpts  = [80] # number of grid points in each segment
#
if gridflag == 0: # medium resolution grid <1E7 entries in 301 table. This is about 250,000 values.
    # define the temperature array for the SESAME table
    # 0 K is changed to 1 K
    # 298 K is inserted into the table if it is not present in the original array
    gridtvals  = [0, 10000., 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [200.,200.,151.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    # 0 g/cm3 is removed if present
    gridrvals  = [1.e-20, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [150, 76, 55, 206, 40] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    # 0 MJ/K/kg is removed if present
    gridsvals  = [0.,0.03] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [1] # 1 for linear, 0 for log in each segment
    gridsnpts  = [500] # number of grid points in each segment
#
if gridflag == 1: # new form for the high resolution table entries in 301 table. Try for 1e6 values
    # define the temperature array for the SESAME table
    # 0 K is changed to 1 K
    # T0REF is inserted into the table if it is not present in the original array
    # need to have extra temperature resolution around the triple point
    gridtvals  = [0, 100, 1600, 2000, 10000., 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,1,1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [2, 151., 140, 300.,100.,50.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    # R0REF is inserted into the table if it is not present in the original array
    # need to have extra density resolution around the triple point
    # 0 g/cm3 is removed if present
    gridrvals  = [1.e-20, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [200, 100, 160, 300, 100] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    # 0 MJ/K/kg is removed if present
    gridsvals  = [1.e-4,.25] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [0] # 1 for linear, 0 for log in each segment
    gridsnpts  = [800] # number of grid points in each segment

if gridflag == 11: # new form for the high resolution table entries in 301 table. Try for 1e6 values
    # define the temperature array for the SESAME table
    # 0 K is changed to 1 K
    # T0REF is inserted into the table if it is not present in the original array
    # need to have extra temperature resolution around the triple point
    gridtvals  = [0, 100, 1600, 2000, 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [5, 201., 200, 200.,50.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    # R0REF is inserted into the table if it is not present in the original array
    # need to have extra density resolution around the triple point
    # 0 g/cm3 is removed if present
    gridrvals  = [1.e-30, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [300, 100, 160, 200, 100] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    # 0 MJ/K/kg is removed if present
    gridsvals  = [1.e-4,.25] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [0] # 1 for linear, 0 for log in each segment
    gridsnpts  = [800] # number of grid points in each segment
#----------------------------------------
# END USER INPUT
#----------------------------------------
In [5]:
# THIS CELL HAS CODE TO GENERATE THE GRIDS AND WRITES TO A FILE TO BE READ BY ANEOS
print('GENERATING ARRAYS FOR THE INDICES OF THE TABULATED EOS.')
#-------------- the following code generates the temperature array
gridtarr = []
for iir in  range(1,len(gridtvals)):
    if gridtstyle[iir-1] == 1:
        # linear region
        temptarr = np.arange(gridtnpts[iir-1])/gridtnpts[iir-1]*(gridtvals[iir]-gridtvals[iir-1])+gridtvals[iir-1]
    if gridtstyle[iir-1] == 0:
        # log region
        logstep = (np.log10(gridtvals[iir])-np.log10(gridtvals[iir-1]))/gridtnpts[iir-1]
        temptarr = np.power(10.,np.log10(gridtvals[iir-1])+(np.arange(gridtnpts[iir-1]))*logstep)
    gridtarr = np.concatenate((gridtarr,temptarr))
# check if T0REF already in the grid
#indexT0REF = get_indexes(T0REF,gridtarr)
indexT0REF = np.where(abs(T0REF-gridtarr)<1.E-4)[0]
if len(indexT0REF) == 0:
    # insert T0REF K into the array
    indexT0REF = np.where((gridtarr-T0REF)>0.)[0]
    gridtarr = np.concatenate((gridtarr[0:indexT0REF[0]],[T0REF],gridtarr[indexT0REF[0]::]))
    if indexT0REF[0]>1:
        print('Inserted the REFERENCE TEMPERATURE INTO THE ARRAY:',gridtarr[indexT0REF[0]-2:indexT0REF[0]+3])
    else:
        print('Inserted the REFERENCE TEMPERATURE INTO THE ARRAY:',gridtarr[indexT0REF[0]-1:indexT0REF[0]+3])
else:
    print('T0REF already in temperature array: ',T0REF,gridrarr[indexT0REF[0]-2:indexT0REF[0]+2])
# Remove 0 K or convert to 1 K
indexzero = get_indexes(0.,gridtarr)
indexone = get_indexes(1.,gridtarr)
if len(indexzero) == 1 and len(indexone) == 0:
    gridtarr[indexzero] = 1.
    print('Zero kelvin changed to 1 kelvin.')
else:
    gridtarr = np.delete(gridtarr, indexzero)
    print('Deleted zero K index.')
#
#gridtsize = len(grdtarr)
print('Number of temperature points: ',len(gridtarr))
#print(gridtarr)
#------------ the following code generates the density array
gridrarr = []
for iir in range(1,len(gridrvals)):
    if gridrstyle[iir-1] == 1:
        # linear region
        temprarr = np.arange(gridrnpts[iir-1])/gridrnpts[iir-1]*(gridrvals[iir]-gridrvals[iir-1])+gridrvals[iir-1]
    if gridrstyle[iir-1] == 0:
        # log region
        logstep = (np.log10(gridrvals[iir])-np.log10(gridrvals[iir-1]))/gridrnpts[iir-1]
        temprarr = np.power(10.,np.log10(gridrvals[iir-1])+(np.arange(gridrnpts[iir-1]))*logstep)
    gridrarr = np.concatenate((gridrarr,temprarr))
# check if R0REF already in the grid
#indexR0REF = get_indexes(R0REF,gridrarr)
indexR0REF = np.where(abs(R0REF-gridrarr)<1.E-4)[0]
if len(indexR0REF) == 0:
    # insert R0REF  into the array
    indexR0REF = np.where((gridrarr-R0REF)>0.)[0]
    gridrarr = np.concatenate((gridrarr[0:indexR0REF[0]],[R0REF],gridrarr[indexR0REF[0]::]))
    print('Inserted the REFERENCE DENSITY INTO THE ARRAY:',gridrarr[indexR0REF[0]-2:indexR0REF[0]+3])
    #print('NOT INSERTING REFERENCE DENSITY INTO THE ARRAY.')
else:
    print('R0REF already in density array: ',R0REF,gridrarr[indexR0REF[0]-2:indexR0REF[0]+2])
# Remove 0 g/cm3 if present
indexzero = get_indexes(0.,gridrarr)
if len(indexzero) == 1:
    gridrarr = np.delete(gridrarr, indexzero)
    print('Deleted zero g/cm3 index.')
#
#gridrsize = len(gridrarr)
print('Number of density points: ',len(gridrarr))
#print(gridrarr)
#------------ the following code generates the specific entropy array for gadget tables
gridsarr = []
for iir in range(1,len(gridsvals)):
    if gridsstyle[iir-1] == 1:
        # linear region
        tempsarr = np.arange(gridsnpts[iir-1])/gridsnpts[iir-1]*(gridsvals[iir]-gridsvals[iir-1])+gridsvals[iir-1]
    if gridsstyle[iir-1] == 0:
        # log region
        logstep = (np.log10(gridsvals[iir])-np.log10(gridsvals[iir-1]))/gridsnpts[iir-1]
        tempsarr = np.power(10.,np.log10(gridsvals[iir-1])+(np.arange(gridsnpts[iir-1]))*logstep)
    gridsarr = np.concatenate((gridsarr,tempsarr))
# Remove 0 MJ/K/kg if present
indexzero = get_indexes(0.,gridsarr)
if len(indexzero) == 1:
    gridsarr[indexzero]=1.e-20
#    gridsarr = np.delete(gridsarr, indexzero)
    print('replaced zero MJ/K/kg index with 1.E-20')
#gridssize = len(gridsarr)
#
print('Number of specific entropy points: ',len(gridsarr))
#print(gridsarr)

print('Number of words in the standard SESAME 301 table (NR, NT, R, T, P, IE, HFE)=',2.+len(gridrarr)+len(gridtarr)+len(gridrarr)*len(gridtarr)*3.)
#-----------------------------------------------------------
#
if DEVMODE:
    print('WRITING tablegrid.txt FOR ANEOS in the local directory.\n')
    # write the table arrays to a file for the ANEOS driver to read in
    tablegridfile = open("tablegrid.txt","w")  
    tablegridfile.write("{:.6e}".format(MATID)+'\n')  
    tablegridfile.write("{:.6e}".format(DATE)+'\n')  
    tablegridfile.write("{:.6e}".format(VERSION)+'\n')  
    tablegridfile.write("{:.6e}".format(FMN)+'\n')  
    tablegridfile.write("{:.6e}".format(FMW)+'\n')  
    tablegridfile.write("{:.6e}".format(R0REF)+'\n')  
    tablegridfile.write("{:.6e}".format(K0REF)+'\n')  
    tablegridfile.write("{:.6e}".format(T0REF)+'\n')  
    tablegridfile.write("{:.6e}".format(len(gridrarr))+'\n')  # number of density points
    tablegridfile.write("{:.6e}".format(len(gridtarr))+'\n')  # number of temperature points
    for i in range(0,len(gridrarr)):
        tablegridfile.write("{:.6e}".format(gridrarr[i])+'\n')
    for i in range(0,len(gridtarr)):
        tablegridfile.write("{:.6e}".format(gridtarr[i])+'\n')
    tablegridfile.close()
#
#-----------------------------------------------------
# Output GADGET2 Table information
if MAKEGADGETTABLE:
    print('Number of words in the extended GADGET2 EOS table (NR, NS, R, S, P, T, IE, CS, HFE, CV, KPA, MDQ)=',2.+len(gridrarr)+len(gridtarr)+len(gridrarr)*len(gridtarr)*8.)
#------------------------------------------------------

#----------------------------------------------------
# Sp. Internal Energy grid for Tillotson
#tillearr = np.power(10.,np.arange(600)/600.*6.-3.) # MJ/kg
#tillearr = np.append(0.,tillearr)
#print('Number of specific internal energy points for Tillotson EOS: ',len(tillearr))
#----------------------------------------------------
# Paranoia check
# make sure all table grid values are unique
if (len(np.unique(gridtarr))!=len(gridtarr)):
    print('NON-UNIQUE VALUES IN TEMP ARRAY')
    stopthecode
if (len(np.unique(gridrarr))!=len(gridrarr)):
    print('NON-UNIQUE VALUES IN DENSITY ARRAY')
    stopthecode
if (len(np.unique(gridsarr))!=len(gridsarr)):
    print('NON-UNIQUE VALUES IN SP. ENTROPY ARRAY')
    stopthecode
#if (len(np.unique(tillearr))!=len(tillearr)):
#    print('NON-UNIQUE VALUES IN SP. ENERGY ARRAY')
#    stopthecode

#----------------------------------------------------

print('\nDone generating EOS grids.\n')
print('\nCHECK THAT THE GRID ARRAYS ARE EVEN IN LENGTH.\n')

# END TABLE GRIDDING INFORMATION
GENERATING ARRAYS FOR THE INDICES OF THE TABULATED EOS.
Inserted the REFERENCE TEMPERATURE INTO THE ARRAY: [286.56716418 294.02985075 298.         301.49253731 308.95522388]
Zero kelvin changed to 1 kelvin.
Number of temperature points:  657
Inserted the REFERENCE DENSITY INTO THE ARRAY: [2.17  2.193 2.2   2.216 2.239]
Number of density points:  861
Number of specific entropy points:  800
Number of words in the standard SESAME 301 table (NR, NT, R, T, P, IE, HFE)= 1698551.0
WRITING tablegrid.txt FOR ANEOS in the local directory.

Number of words in the extended GADGET2 EOS table (NR, NS, R, S, P, T, IE, CS, HFE, CV, KPA, MDQ)= 4526936.0

Done generating EOS grids.


CHECK THAT THE GRID ARRAYS ARE EVEN IN LENGTH.

In [6]:
if DEVMODE:
    # call aneos in the local directory
    # output file is in the local directory
    print('RUNNING IN DEVELOPMENT MODE')
    print('Calling aneos in the local directory.')
    completed = subprocess.run(['./aneos-mbam3'])
    print('Aneos call complete. CHECK THAT A NEW ANEOS.OUTPUT FILE WAS CREATED.')
    print('aneos call result: ', completed.returncode)
    print('Check ANEOS.OUTPUT file for errors.')
else:
    print('RUNNING IN RELEASE MODE: NO ANEOS SYSTEM CALL.')
    print('READING IN EOS DISTRIBUTION FILES.')
RUNNING IN DEVELOPMENT MODE
Calling aneos in the local directory.
Aneos call complete. CHECK THAT A NEW ANEOS.OUTPUT FILE WAS CREATED.
aneos call result:  0
Check ANEOS.OUTPUT file for errors.
In [7]:
# READ IN NEW ANEOS MODEL and fill the extEOStable class object
# source in eostable.py
#------------------------------------------------------------------
NewEOS  = etab.extEOStable() # FIRST make new empty EOS object
NewEOS.loadextsesame('NEW-SESAME-EXT.TXT') # LOAD THE EXTENDED 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
NewEOS.loadstdsesame('NEW-SESAME-STD.TXT') # LOAD THE STANDARD 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
NewEOS.MODELNAME = MODELNAME # string set above in user input
NewEOS.MDQ = np.zeros((NewEOS.NT,NewEOS.ND)) # makes the empty MDQ array
#print(NewEOS.units) # these are the default units for SESAME rho-T tables
#'Units: g/cm3, K, GPa, MJ/kg, MJ/kg, MJ/K/kg, cm/s, MJ/K/kg, KPA flag. 2D arrays are (NT,ND).'

# Add the header info to the table. This could be done during the loading. 
# if made from this notebook, these values are set in the user-input above.
# ** MAKE SURE THEY MATCH ANEOS.INPUT **
NewEOS.MATID   = MATID
NewEOS.DATE    = DATE
NewEOS.VERSION = VERSION
NewEOS.FMN     = FMN
NewEOS.FMW     = FMW
NewEOS.R0REF   = R0REF
NewEOS.K0REF   = K0REF
NewEOS.T0REF   = T0REF
NewEOS.P0REF   = P0REF
#
# Load the information from ANEOS.INPUT and ANEOS.OUTPUT
NewEOS.loadaneos(aneosinfname='ANEOS.INPUT',aneosoutfname='ANEOS.OUTPUT')
#
# calculate the principal Hugoniot from the table
if DEVMODE:
    # calculate the Hugoniot and write to a file. This is from the ANEOS REFERENCE STATE
    NewEOS.calchugoniot(r0=NewEOS.R0REF,t0=NewEOS.T0REF,writefilename='NEW-SESAME-HUG.TXT')
    # calculate the Hugoniot and write to a file pyrolitic Glass initial state, which is not on this EOS surface (proxy crystalline solid + liquid + vapor)
    # this result is basically the same as the requested Hugoniot at r0 and t0, where ANEOS interprets the initial density to be in tension. Differences at very low pressures.
    #NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.1565,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
    NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
else:
    # just calculate the Hugoniot and put it in the EOS object
    NewEOS.calchugoniot(r0=NewEOS.R0REF,t0=NewEOS.T0REF)
#    NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.1565,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
    NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
#
# 9390/.06/1.e6=0.1565 Richet 9390 J/mol -> MJ/kg
# calculate the 1-bar profile; loop over temp
NewEOS.onebar.T = np.zeros(NewEOS.NT)
NewEOS.onebar.S = np.zeros(NewEOS.NT)
NewEOS.onebar.rho = np.zeros(NewEOS.NT)
it0 = np.where(NewEOS.T >= NewEOS.T0REF)[0]
id0 = np.arange(NewEOS.ND)#np.where(NewEOS.rho >= 0.8*NewEOS.R0REF)[0]
for iit in range(0,NewEOS.NT):
    NewEOS.onebar.T[iit] = NewEOS.T[iit]
    NewEOS.onebar.S[iit] = np.interp(1.E-4,NewEOS.P[iit,id0],NewEOS.S[iit,id0])
    NewEOS.onebar.rho[iit] = np.interp(1.E-4,NewEOS.P[iit,id0],NewEOS.rho[id0])
    #print(iit,NewEOS.onebar.rho[iit],NewEOS.onebar.T[iit],NewEOS.onebar.S[iit]*1.E3)
Done loading ANEOS files.
ANEOS WAS CALLED WITH THE FOLLOWING INPUT, LOADED FROM FILE  ANEOS.INPUT
 *         EOS#      string for EOS name           .....Hugoniot initial........ long
 *                                                 T(eV)           Rho(g/cc)     flag
 * Reference solid Hugoniot at STP initial density and temperature
 ANEOS1    -1        'FUSED_SILICA_v0.2'           THUG=-1         RHUG=2.2      LONG

 *
 *         nelem     type      rho0      temp0(K)  press0    B0/-Cbulk Gamma0    Tdebye(K)
 *         V01       V02       V03       V04       V05       V06       V07       V08
 ANEOS2    2         4         2.20      298.      1.e6      0.40E12   0.57      -1600.

rho0= 2.2
 *         TG model  Gamma     Esep      Tmelt(K)  ..Critical Pt..     Thermal Conduction
 *         or S1     3*C24     erg/g     or -Em    C53       C54       H0        C41
 *         V09       V10       V11       V12       V13       V14       V15       V16
 ANEOS3    1.        6.0       1.20E11   1996.     -10.0E11  0.90      0.        0.

 *         solid     .....Solid-solid phase transition data.....                 rhol/rhos
 *         rhomin    D1        D2        D3        D4        D5        Hfusion   or -rholiq
 *         V17       V18       V19       V20       V21       V22       V23       V24
 ANEOS4    0.        2.92      4.30      0.4E11    2.85E12   1.E1      0.66E10   0.93

 *         Cold comp. limits    ...Liquid EOS params...       ...Interpolation params...
 *         Up        L0        alpha     beta      gamma     C60       C61       C62
 *         V25       V26       V27       V28       V29       V30       V31       V32
 ANEOS5    0.        0.        0.        0.        0.        0.0       -0.851     0.5

 *         Ioniz.    Reactive_chemistry  ...............Molecular Clusters...............
 *         model     Eshift    Sshift    Atoms     Ebind(eV) RotDOF    Rbond     VibDOF
 *         V33       V34       V35       V36       V37       V38       V39       V40
 ANEOS6    0         0.        0.        2         5.0       2.0       1.5E-8    1

 *         Theta_vib Pcold/Expanded      fcv*3nR   QCC1      QCC6      Unused    Unused
 *         (K)       Mieflag   a_exp     adjust    rhomin    psitrans  
 *         V41       V42       V43       V44       V45       V46       V47       V48
 ANEOS7    2000.     1         1.70      1.29      1.E-30    1.D5      0.        0.

 *
 *         Elements in material (number must = nelem)
 *         Z         Atomic Fraction
 ANEOS8    8         0.66667

 ANEOS10   14        0.33333

READING DATA FROM ANEOS OUTPUT FILE  ANEOS.OUTPUT
Writing Hugoniot to file:  NEW-SESAME-HUG.TXT
r1= 2.65
468 32 2.2 298.0 0.0001 0.0 0.0 1.4440455300000001
Writing Hugoniot to file:  NEW-SESAME-HUG-GLASS.TXT

May 2023 No tension table using adjacent values. See no tension notebook for a little bit of actual interpolation.

In [8]:
# Make a NOTENSION TABLE and GADGET TABLE
################### REMOVE TENSION REGION ###########################
## The GADGET2 tables do not have a tension region. 
## If the gadget2 table is selected, then must generate a no-tension table
if MAKEGADGETTABLE:
    REMOVETENSIONFLAG=1  # 0- keep tension region; 1- remove tension region

# minimum pressure 1.E-30
# minimum sp energy about 1.E-4 before negative entropy region
# minimum sp entropy about 1.E-5 before negative entropy region
# minimum HFE about 1.E-3 before negative entropy region
# minimum sound speed about 1.E-10 before negative entropy region
# minimum sp heat capacity about 1.E-4 before negative entropy region

## FUTURE WORK: THIS SECTION COULD USE THE LEVER RULE TO FILL IN THE MISSING POINTS
## SINCE MANY VALUES ARE NEAR THE LOWER LIMITS IN THE SUBLIMATION REGION ANYWAY, 
## HERE, TENSION REGION VALUES ARE FILLED IN WITH THE ADJACENT S-V DOME VALUES
## This is obviously bad for interpolations. Fix for future work.
if REMOVETENSIONFLAG:
    # First make a deep copy of the original table with the tension region
    NewEOSorig = deepcopy(NewEOS)
    NewEOS.MODELNAME = NewEOS.MODELNAME+'notension'
    # THIS SECTION REWRITES NEGATIVE PRESSURE REGION
    # The tension region is in the solid-vapor phase boundary. 
    # Use the liquid-vapor phase boundary. Below the temperature of the triple point, it is the sublimation curve.
    # But the pressures for the solid calculation is spurious at very low pressures, so use the vapor pressure.
    # HFE = GFE + TS
    print('REMOVING TENSION REGION FROM THE ORIGINAL ANEOS TABLE.')
    for itt in range(NewEOS.NT-1,-1,-1):
        tmp = np.where(NewEOS.P[itt,:] < 0.)
        if (len(tmp[0]) > 0):
            badind = tmp[0]
            if badind[0] > 0:
                #print(itt, len(badind), NewEOS.T[itt], 'tension only', NewEOS.P[itt,badind[0]-4:badind[0]+1])
                # at this temperature, all densities are the same pressure in the vapor dome
                # so use the adjacent positive pressure value
                # do the same for all the variables since there is very little change in this region of the dome.
                # Errors introduced by interpolating the vapor curve are comparable to
                # recalculating the lever rule. Would need to call ANEOS again to get good end point values
                # to calculate the mass-weighted value via the lever rule.
                NewEOS.P[itt,badind] = np.full(len(badind),NewEOS.P[itt,badind[0]-1])
                NewEOS.U[itt,badind] = np.full(len(badind),NewEOS.U[itt,badind[0]-1])
                NewEOS.S[itt,badind] = np.full(len(badind),NewEOS.S[itt,badind[0]-1])
                NewEOS.A[itt,badind] = np.full(len(badind),NewEOS.A[itt,badind[0]-1])
                NewEOS.cs[itt,badind] = np.full(len(badind),NewEOS.cs[itt,badind[0]-1])
                NewEOS.cv[itt,badind] = np.full(len(badind),NewEOS.cv[itt,badind[0]-1])
                NewEOS.KPA[itt,badind] = np.full(len(badind),2)
                NewEOS.MDQ[itt,badind] = np.full(len(badind),-1) # replaced tension region flag
            else:
                print('SOMETHING VERY WRONG WITH TENSION IN THE TABLE')
                makethecodestop # please tell me a more graceful way to stop the notebook
    #
    # THIS SECTION RESETS NEGATIVE SPECIFIC ENTROPIES WITH ADJACENT VALUES OR MINIMUM POSITIVE VALUES
    # THESE VALUES TYPICALLY DEFINE THE COLD COMPRESSION LIMIT OF THE MATERIAL
    # OR SPURIOUS NEGATIVE ENTROPIES AT VERY LOW PRESSURES (IMPERFECT DEBYE MODEL IN ANEOS)
    for itt in range(NewEOS.NT-1,-1,-1):
        tmp = np.where(NewEOS.S[itt,:] < 0.)
        if (len(tmp[0]) > 0):
            badind = tmp[0]
            if badind[0] > 0:
                #print(itt, antarr[itt], andarr[badind[0]], badind[0],len(badind),'negative entropy', ansarr[itt,badind[0]],ansarr[itt,badind[0]-1])
                NewEOS.P[itt,badind] = np.full(len(badind),NewEOS.P[itt,badind[0]-1]) # use the adjacent pressure at this temperature in the dome
                NewEOS.U[itt,badind] = np.full(len(badind),NewEOS.U[itt,badind[0]-1])
                NewEOS.S[itt,badind] = np.full(len(badind),NewEOS.S[itt,badind[0]-1])
                NewEOS.A[itt,badind] = np.full(len(badind),NewEOS.A[itt,badind[0]-1])
                NewEOS.cs[itt,badind] = np.full(len(badind),NewEOS.cs[itt,badind[0]-1])
                NewEOS.cv[itt,badind] = np.full(len(badind),NewEOS.cv[itt,badind[0]-1])
                #NewEOS.KPA[itt,badind]  # leave as original region flag
                NewEOS.MDQ[itt,badind] = np.full(len(badind),-2) # negative entropies flag
                #print('Some NEGATIVE  S AT THIS INDEX, TEMP(K): ',itt, NewEOS.T[itt],NewEOS.P[itt,badind[0]-1],NewEOS.P[itt,badind[0]-2]) #, badind[0],len(badind),'negative entropy', NewEOS.S[itt,badind[0]])
            if badind[0] == 0:
                print('All S VALUES NEGATIVE AT THIS INDEX, TEMP(K): ',itt, NewEOS.T[itt]) #, badind[0],len(badind),'negative entropy', NewEOS.S[itt,badind[0]])
                NewEOS.S[itt,badind] = np.full(len(badind),1.E-10)            
                NewEOS.U[itt,badind] = np.full(len(badind),1.E-10) # also replace the negative sp. internal energies           
                NewEOS.MDQ[itt,badind] = np.full(len(badind),-2) # negative entropies flag

    #### ONLY WRITE A NEW TABLE IN DEVELOPMENT MODE
    if DEVMODE:
        #### WRITE NEW SESAME TABLE WITH NO TENSION REGION
        #     WRITE STANDARD (SHORT) SESAME FILE
        #     WRITE SESAME TABLE TO FILE LIMITED TO P, E, HFE
        NewEOS.writestdsesame(writestdsesfname="NEW-SESAME-STD-NOTENSION.TXT")
        NewEOS.writeextsesame(writeextsesfname="NEW-SESAME-EXT-NOTENSION.TXT")
    if DEVMODE:
        #### WRITE NEW SESAME TABLE WITH NO TENSION REGION
        #     WRITE STANDARD (SHORT) SESAME FILE
        #     WRITE SESAME TABLE TO FILE LIMITED TO P, E, S
        NewEOS.writegadgetinitsesame(writegadgetinitsesfname="NEW-SESAME-GADGETINIT-NOTENSION.TXT")

#### ONLY WRITE A NEW TABLE IN DEVELOPMENT MODE
if DEVMODE:
    #### WRITE MDQ Flag as a 301-style table
    NewEOS.writemdqsesame(writemdqsesfname="NEW-SESAME-MDQ.TXT")
#
#--------------------------------------------------------------------------
########## MAKE GADGET2 EOS TABLE ARRAYS ###################
# LATER WE CAN TURN THIS INTO A FUNCTION IN THE GADGET EOS OBJECT
# FLAG IS SET ABOVE IN THE USER INPUT SECTION
# MAKEGADGETTABLE = 1 # 1 - do it ; 0 - skip it
if MAKEGADGETTABLE:
    print('Starting interpolation for GADGET table.')
    GadEOS  = etab.extGADtable() # make new empty Gadget EOS object
    #
    # rho array is andarr to match the SESAME table
    # S array is gridsarr
    GadEOS.ND  = len(gridrarr)
    GadEOS.NS  = len(gridsarr)
    GadEOS.rho = gridrarr
    GadEOS.S   = gridsarr
    GadEOS.P   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.T   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.U   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.A   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.cs  = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.cv  = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.KPA = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.MDQ = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.MODELNAME = NewEOS.MODELNAME+'-gadget-linear-int'
    #
    # some combinations of density and sp. entropy are not in the EOS surface
    # in these cases, apply minimum values
    gadpmin  = np.min(NewEOS.P) # GPa
    gadsmin  = 1.E-30     # MJ/K/kg
    gademin  = np.min(NewEOS.U)    # MJ/kg
    gadtmin  = 1.     # K
    gadcmin  = np.min(NewEOS.cs) # cm/s
    gadcvmin = 0.     # MJ/K/kg
    gadamin  = np.min(NewEOS.A) # can be negative    # MJ/kg
    print('Minimum values used in the rational function interpolation scheme:')
    print('gadpmin,gadsmin,gademin,gadtmin,gadcmin,gadcvmin,gadamin=')
    print(gadpmin,gadsmin,gademin,gadtmin,gadcmin,gadcvmin,gadamin)
    #
    # some combinations of density and sp. entropy are not in the EOS surface
    # in these cases, fill the region with the adjacent real values from the original T-RHO table
    # ----------------------------------------------------------
    # SIMPLE LINEAR INTERPOLATION TO FLIP T-RHO to S-RHO table
    print("LINEAR INTERPOLATION FOR GADGET TABLE.")
    for iid in range(0,GadEOS.ND):
        # same density arrays for the T-rho table and S-rho table
        sall1d = NewEOS.S[:,iid]
        pall1d = NewEOS.P[:,iid]
        eall1d = NewEOS.U[:,iid]
        tall1d = NewEOS.T[:]
        call1d = NewEOS.cs[:,iid]
        cvall1d = NewEOS.cv[:,iid]
        aall1d = NewEOS.A[:,iid]
        kall1d = NewEOS.KPA[:,iid]
        mdqall1d = NewEOS.MDQ[:,iid]
        # interpolate P at the new S
        testinterp = interpolate.interp1d(sall1d,pall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.P[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.P[tmpind,iid] = pall1d[NewEOS.NT-1]
        # interpolate E at the new S
        testinterp = interpolate.interp1d(sall1d,eall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.U[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.U[tmpind,iid] = eall1d[NewEOS.NT-1]
        # interpolate T at the new S
        testinterp = interpolate.interp1d(sall1d,tall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.T[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.T[tmpind,iid] = tall1d[NewEOS.NT-1]
        # interpolate sound speed at the new S
        testinterp = interpolate.interp1d(sall1d,call1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.cs[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.cs[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate sp. heat capacity at the new S
        testinterp = interpolate.interp1d(sall1d,cvall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.cv[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.cv[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate HFE at the new S
        testinterp = interpolate.interp1d(sall1d,aall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.A[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.A[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate phase information KPA flag at the new S
        testinterp = interpolate.interp1d(sall1d,kall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.KPA[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.KPA[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate model development quality flag at the new S
        testinterp = interpolate.interp1d(sall1d,mdqall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.MDQ[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.MDQ[tmpind,iid] = call1d[NewEOS.NT-1]
    #
    GadEOSlin = deepcopy(GadEOS)
    #
    '''
    # ----------------------------------------------------------
    # 1D RATIONAL FUNCTION INTERPOLATION METHOD FROM KERLEY
    print("RATIONAL FUNCTION INTERPOLATION TO GENERATE GADGET2 EOS TABLE.")
    for iid in range(0,GadEOS.ND):
        GadEOS.MODELNAME = NewEOS.MODELNAME+'-gadget-ratfunc-int'
        #
        # same density arrays for the T-rho table and S-rho table
        sall1d  = NewEOS.S[:,iid]
        pall1d  = NewEOS.P[:,iid]
        eall1d  = NewEOS.U[:,iid]
        tall1d  = NewEOS.T[:]
        call1d  = NewEOS.cs[:,iid]
        cvall1d = NewEOS.cv[:,iid]
        aall1d  = NewEOS.A[:,iid]
        #kall1d = NewEOS.KPA[:,iid]
        #mdqall1d = NewEOS.MDQ[:,iid]
        #
        # interpolate P at the new S
        GadEOS.P[:,iid] = rf1d(sall1d,pall1d,gridsarr,gadpmin)
        # interpolate E at the new S
        GadEOS.U[:,iid] = rf1d(sall1d,eall1d,gridsarr,gademin)
        # interpolate T at the new S
        GadEOS.T[:,iid] = rf1d(sall1d,tall1d,gridsarr,gadtmin)
        # interpolate sound speed at the new S
        GadEOS.cs[:,iid] = rf1d(sall1d,call1d,gridsarr,gadcmin)
        # interpolate sp. heat capacity at the new S
        GadEOS.cv[:,iid] = rf1d(sall1d,cvall1d,gridsarr,gadcvmin)
        # interpolate HFE at the new S
        GadEOS.A[:,iid] = rf1d(sall1d,aall1d,gridsarr,gadamin)
        # Do not interpolate KPA or MDQ with rational function interp. Use Linear.
        # do not overwrite the linear interpolation above.
    '''
    #
    print('Done interpolating the gadget2 table.')
    #print('Using RATIONAL FUNCTION INTERPOLATION EXCEPT FOR KPA and MDQ FLAGS (bilinear).')
    if DEVMODE:
        # ONLY WRITE GADGET2 TABLE TO DISK IN DEVELOPMENT MODE
        #### WRITE NEW GADGET2 TABLE WITH NO TENSION REGION
        #### This ascii table format follows the sesame table format of 5 columns of numbers in 16.8e format
        #### the STYLE variable is just a counter to 5 to know when to write a newline
        GadEOS.writestdgadget(writestdgadgetfname="NEW-GADGET-STD-NOTENSION.TXT")
        GadEOS.writeextgadget(writeextgadgetfname="NEW-GADGET-EXT-NOTENSION.TXT")
REMOVING TENSION REGION FROM THE ORIGINAL ANEOS TABLE.
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  10 137.3134
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  9 129.8507
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  8 122.3881
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  7 114.9254
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  6 107.4627
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  5 100.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  4 80.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  3 60.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  2 40.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  1 20.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  0 1.0
Done writing the STD SESAME 301 table to local directory:  NEW-SESAME-STD-NOTENSION.TXT
Done writing the EXT SESAME 301 table to local directory:  NEW-SESAME-EXT-NOTENSION.TXT
Done writing the GADGET INIT SESAME 301 notension table to local directory:  NEW-SESAME-GADGETINIT-NOTENSION.TXT
Done writing the MDQ Flag as a 301-style table to local directory:  NEW-SESAME-MDQ.TXT
Starting interpolation for GADGET table.
Minimum values used in the rational function interpolation scheme:
gadpmin,gadsmin,gademin,gadtmin,gadcmin,gadcvmin,gadamin=
9.99999968e-31 1e-30 1e-10 1.0 1e-10 0.0 -3116194.75
LINEAR INTERPOLATION FOR GADGET TABLE.
Done interpolating the gadget2 table.
Done writing the STD GADGET TABLE FILE with notension to local directory:  NEW-GADGET-STD-NOTENSION.TXT
Done writing the EXTENDED GADGET TABLE FILE with notension to local directory:  NEW-GADGET-EXT-NOTENSION.TXT

Read in silica data¶

In [9]:
%%capture 
## This cell loads the modules, data, and EOS parameters
## capture magic function hides the output; comment out to see diagnostic info during data load

# These scripts and data will be distributed in a separate GitHub repository (Amodeo et al., in prep.)

qtz_liq_hug = 'KD09' # or `KD13` # pick which quartz Hugoniot equation
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_paths
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_modules
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_silica_data
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/define_eos
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/define_mc
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_vapor_domes.py
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_liquid_model.py
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_solid_models.py
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/hugoniot-calcs.py
In [10]:
# Read in additional SiO2 system data that is not read in load_silica_data.py
# silicadatadir = SET IN LOAD_SILICA_DATA.PY

# Use these units:
units_str = 'rho g/cm3, T K, P GPa, U MJ/kg, A MJ/kg, S kJ/K/kg, cs cm/s, cv MJ/K/kg, us km/s, up km/s'

# Bourova, E., Richet, P., 1998. Quartz and Cristobalite: high-temperature cell parameters and volumes of fusion. Geophys. Res. Lett. 25, 2333–2336. https://doi.org/10.1029/98GL01581
# Bourova & Richet 1998 volumes of fusion qtz, cristobalite, liquid
qtz_mc_BR98_file = silicadatadir+'BourovaRichet1998-qtz-melt-curve.txt'
qtz_mc_BR98 = etab.EOSmeltcurve()
qtz_mc_BR98.Pl = np.loadtxt(qtz_mc_BR98_file,delimiter=',',skiprows=2,usecols=[0])*1.E8/1.e9 # kbar -> GPa
qtz_mc_BR98.Ps = np.loadtxt(qtz_mc_BR98_file,delimiter=',',skiprows=2,usecols=[0])*1.E8/1.e9 # kbar -> GPa
qtz_mc_BR98.T = np.loadtxt(qtz_mc_BR98_file,delimiter=',',skiprows=2,usecols=[1]) # K
qtz_mc_BR98.units = units_str
qtz_mc_BR98.label = 'Melt curve Bourova & Richet 1998'
#print(qtz_mc_BR98.Pl)

# Millot, M., Dubrovinskaia, N., Černok, A., Blaha, S., Dubrovinsky, L., Braun, D.G., Celliers, P.M., Collins, G.W., Eggert, J.H., Jeanloz, R., 2015. Shock compression of stishovite and melting of silica at planetary interior conditions. Science (80-. ). 347, 418–420. https://doi.org/10.1126/science.1261507
# Millot et al. Science 2015 Melt curve fitted Simon's equation using stishovite crossing point
#stv_mc_MM15_file   = silicadatadir+'Millot-melt-curve.txt'
#stv_mc_MM15 = EOSmeltcurve()
#stv_mc_MM15.Pl = np.loadtxt(stvmcfile,delimiter=',',skiprows=0,usecols=[0])*1.e9 # Pa
#stv_mc_MM15.Ps = stv_mc.Pl
#stv_mc_MM15.T = np.loadtxt(stvmcfile,delimiter=',',skiprows=0,usecols=[1])*1.e3 # K
stv_mc_MM15 = etab.EOSmeltcurve()
stv_mc_MM15.Pl = np.arange(0,600) # GPa
stv_mc_MM15.Ps = np.arange(0,600) # GPa
stv_mc_MM15.T = 1968.5 + 307.8*np.power(stv_mc_MM15.Pl,0.485) # K with P in GPa
stv_mc_MM15.units = units_str
stv_mc_MM15.label = 'Melt curve Millot et al. 2015'
#print(stv_mc_MM15.T)

# McCoy, C.A., Gregor, M.C., Polsin, D.N., Fratanduono, D.E., Celliers, P.M., Boehly, T.R., Meyerhofer, D.D., 2016. Measurements of the sound velocity of shock-compressed liquid silica to 1100 GPa. J. Appl. Phys. 120. https://doi.org/10.1063/1.4972338
# McCoy et al. JAP 2016 Gammas from sound speeds on the Hugoniot
qtz_hug_McCoy16_file = silicadatadir+'McCoy2016-qtz-hug-gamma.txt'
qtz_hug_McCoy16 = etab.EOShugoniot()
qtz_hug_McCoy16.P = np.loadtxt(qtz_hug_McCoy16_file,delimiter=' ',skiprows=2,usecols=[1]) # GPa
qtz_hug_McCoy16.cs = np.loadtxt(qtz_hug_McCoy16_file,delimiter=' ',skiprows=2,usecols=[5]) # km/s
qtz_hug_McCoy16.gamma = np.loadtxt(qtz_hug_McCoy16_file,delimiter=' ',skiprows=1,usecols=[7]) # [-] 
qtz_hug_McCoy16.rho = np.interp(qtz_hug_McCoy16.P,qtz_hug_KD.P,qtz_hug_KD.rho) # g/cm3
qtz_hug_McCoy16.units = units_str
qtz_hug_McCoy16.label = 'Qtz McCoy et al. 2016'
#print(mccoy_hug.S)
#print(mccoy_hug.rho)

# NIST-JANAF web site SiO2
qtz_1bar_JANAF_file = silicadatadir+'JANAF-qtz-crystal-liquid.txt'
qtz_1bar_JANAF = etab.EOS1barcurve()
qtz_1bar_JANAF.T = np.loadtxt(qtz_1bar_JANAF_file,delimiter='\t',skiprows=2,usecols=[0]) # K
qtz_1bar_JANAF.S = np.loadtxt(qtz_1bar_JANAF_file,delimiter='\t',skiprows=2,usecols=[2])/0.0600843/1.e3 # J/K/mol -> kJ/K/kg
qtz_1bar_JANAF.units = units_str
qtz_1bar_JANAF.label = '1 bar JANAF'
#print(qtz_1bar_JANAF.S)
#janafgasfile = silicadatadir+'janaf-gas.txt'
#janafgas_1bar = EOS1barcurve()
#janafgas_1bar.T = np.loadtxt(janafgasfile,delimiter='\t',skiprows=3,usecols=[0]) # K
#janafgas_1bar.S = np.loadtxt(janafgasfile,delimiter='\t',skiprows=3,usecols=[2])/0.0600843 # J/K/mol -> J/K/kg
#print(janafgas_1bar.S) # this does not seem to be right?? does not agree with Melosh 2007

## read in Connolly data
C16highs = np.loadtxt(silicadatadir+'Connolly-high-cp.csv',delimiter=',',usecols=[0])
C16hight = np.loadtxt(silicadatadir+'Connolly-high-cp.csv',delimiter=',',usecols=[1])
C16bsrho  = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-T-rho.csv',delimiter=',',usecols=[0])/1.e3
C16bstemp = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-T-rho.csv',delimiter=',',usecols=[1])


C16lows = np.loadtxt(silicadatadir+'Connolly-low-cp.csv',delimiter=',',usecols=[0])
C16lowt = np.loadtxt(silicadatadir+'Connolly-low-cp.csv',delimiter=',',usecols=[1])
C16bdrho  = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-dashed-T-rho.csv',delimiter=',',usecols=[0])/1.e3
C16bdtemp = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-dashed-T-rho.csv',delimiter=',',usecols=[1])

C16sils = np.loadtxt(silicadatadir+'Connolly-SiL-VWRK.csv',delimiter=',',usecols=[0])
C16silt = np.loadtxt(silicadatadir+'Connolly-SiL-VWRK.csv',delimiter=',',usecols=[1])
C16rsrho  = np.loadtxt(silicadatadir+'Connolly-Fig4b-red-T-rho.csv',delimiter=',',usecols=[0])/1.e3
C16rstemp = np.loadtxt(silicadatadir+'Connolly-Fig4b-red-T-rho.csv',delimiter=',',usecols=[1])

# Load Karki gammas 2007
v0 = 45.8 # A/SiO2
rho0 = sio2mol/v0*1.6605391 # g/cm3

# STS estimated melt curve line for data review draw as dotted above 155 GPa
dat_file   = silicadatadir+'karki-gamma.txt'
dat         = etab.EOSmeltcurve()
dat.matid    = [174]
dat.material = ['Quartz']
dat.formula  = ['SiO2']
dat.rho_ref = [2.650] # g/cm3
sio2mol = 28.085+ 2*15.999
vol=np.loadtxt(dat_file,delimiter=',',skiprows=1,usecols=[0])*v0
dat.rho      = sio2mol/vol*1.6605391 # amu/a^3 -> g/cm3
dat.rho_units = ['g/cm3']
dat.gamma=np.loadtxt(dat_file,delimiter=',',skiprows=1,usecols=[1])

dat_file   = silicadatadir+'karki-cv.txt'
dat.cv=np.loadtxt(dat_file,delimiter=',',skiprows=1,usecols=[1])*3*8.314/(sio2mol/1000) # J/K/kg

#print(dat.gamma)
#print(dat.rho)
#print(dat.cv)
karkigcvorig = deepcopy(dat)
dat=[]
vol=[]

Color mesh plots of the SESAME Rho-T table¶

In [11]:
# PLOT COLOR SHADED VALUES RHO-T TABLE
lw=1
sfont = 15
ds = (10,10)
plt.rcParams['pcolor.shading'] ='nearest'
In [12]:
%matplotlib inline
fig, ([ax0, ax1], [ax2,ax3], [ax4,ax5])= plt.subplots(nrows=3,ncols=2,figsize=(14,24))
plt.subplots_adjust(wspace=0.4)
#plt.tight_layout(pad=0.6)

##================ Pressure
im1=ax0.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.P)+25)/29.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax0.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax0.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax0.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax0.set_ylim(0.,20000)
ax0.set_xlim(0,30.)
ax0.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax0.set_ylabel('Temperature (K)',size=sfont)
ax0.set_title('Colors=log$_{10}$(P (GPa))',size=sfont)
ax0.tick_params(labelsize=sfont)

divider0 = make_axes_locatable(ax0)
cax0 = divider0.append_axes('right', size='5%', pad=0.05)
cbar0 = fig.colorbar(im1, cax=cax0, ticks = [],orientation='vertical')
#cbar0.ax.set_yticklabels(['  ','  ','  '])  # vertical colorbar
cbar0.set_label('Log$_{10}$(P) [GPa]',size=sfont)
#cbar0.ax.set_yticklabels(['-25','-11','+4'])  # vertical colorbar

im1=ax1.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.P)+25)/29.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax1.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax1.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax1.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax1.set_ylim(10.,1.E7)
ax1.set_xlim(1.E-21,100.)
ax1.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax1.set_ylabel('Temperature (K)',size=sfont)
ax1.set_title('Colors=log$_{10}$(P (GPa))',size=sfont)
ax1.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-25','-11', '+4'])  # vertical colorbar
cbar.set_label('Log$_{10}$(P) [GPa]',size=sfont)


#================ Sp. Int. Energy
im1=ax2.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.U)+4.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax2.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax2.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax2.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax2.set_ylim(0.,20000)
ax2.set_xlim(0,30.)
ax2.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax2.set_ylabel('Temperature (K)',size=sfont)
ax2.set_title('Colors=log$_{10}$(Sp. Int. Energy (MJ/kg))',size=sfont)
ax2.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
# does not want to label ticks in this location
cbar = fig.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['','',''])  # vertical colorbar
cbar.set_label('Log$_{10}$(E) [MJ/kg]',size=sfont)
#cbar.ax.set_yticklabels(['-4','+1','+6'])  # vertical colorbar


ax3.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.U)+4.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax3.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax3.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax3.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax3.set_ylim(10.,1.E7)
ax3.set_xlim(1.E-21,100.)
ax3.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax3.set_ylabel('Temperature (K)',size=sfont)
ax3.set_title('Colors=log$_{10}$(E (MJ/kg))',size=sfont)
ax3.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax3)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-4', '+1','+6'])  # vertical colorbar
cbar.set_label('Log$_{10}$(E) [MJ/kg]',size=sfont)

#================ Sp. Entropy
im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.S)+4.)/4.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax4.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax4.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax4.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax4.set_ylim(0.,20000)
ax4.set_xlim(0,30.)
ax4.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax4.set_ylabel('Temperature (K)',size=sfont)
ax4.set_title('Colors=log$_{10}$(Sp. Entropy (MJ/K/kg))',size=sfont)
ax4.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax4)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['-4', '-2', '0'])  # vertical colorbar
cbar.set_label('Log$_{10}$(S) [MJ/K/kg]',size=sfont)


ax5.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.S)+4.)/4.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax5.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax5.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax5.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax5.set_ylim(10.,1.E7)
ax5.set_xlim(1.E-21,100.)
ax5.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax5.set_ylabel('Temperature (K)',size=sfont)
ax5.set_title('Colors=log$_{10}$(Sp. Entropy (MJ/K/kg))',size=sfont)
ax5.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax5)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-4', '-2','0'])  # vertical colorbar
cbar.set_label('Log$_{10}$(S) [MJ/K/kg]',size=sfont)


##============= start new plot
fig2, ([ax0, ax1], [ax2,ax3], [ax4,ax5])= plt.subplots(nrows=3,ncols=2,figsize=(14,24))
plt.subplots_adjust(wspace=0.4)

#================ Sound Speed
im1=ax0.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.cs)+3.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax0.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax0.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax0.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax0.set_ylim(0.,20000)
ax0.set_xlim(0,30.)
ax0.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax0.set_ylabel('Temperature (K)',size=sfont)
ax0.set_title('Colors=log$_{10}$(Sound Speed (cm/s))',size=sfont)
ax0.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax0)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['-3', '+2', '+7'])  # vertical colorbar
cbar.set_label('Log$_{10}$(Cs) [cm/s]',size=sfont)


ax1.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.cs)+3.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax1.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax1.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax1.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax1.set_ylim(10.,1.E7)
ax1.set_xlim(1.E-21,100.)
ax1.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax1.set_ylabel('Temperature (K)',size=sfont)
ax1.set_title('Colors=log$_{10}$(sound speed (cm/s))',size=sfont)
ax1.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig2.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-3', '+2', '+7'])  # vertical colorbar
cbar.set_label('Log$_{10}$(Cs) [cm/s]',size=sfont)


#================ HFE
amin = np.min(NewEOS.A)
amin = -1000.
alogmax = 1500. #np.max(anaarrodd)-amin #np.max(np.log10(anaarrodd-amin))
#print('amin, log10(amin)',amin,np.log10(-amin))
#print('alogmax=',alogmax)
im1=ax2.pcolormesh(NewEOS.rho,NewEOS.T,(NewEOS.A-amin)/alogmax,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax2.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax2.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax2.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax2.set_ylim(0.,20000)
ax2.set_xlim(0,30.)
ax2.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax2.set_ylabel('Temperature (K)',size=sfont)
ax2.set_title('Colors=HFE (MJ/kg)',size=sfont)
ax2.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['-1000','-500','0','500'])  # vertical colorbar
cbar.set_label('Log$_{10}$(HFE) [MJ/kg]',size=sfont)


ax3.pcolormesh(NewEOS.rho,NewEOS.T,(NewEOS.A-amin)/alogmax,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax3.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax3.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax3.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax3.set_ylim(10.,1.E7)
ax3.set_xlim(1.E-21,100.)
ax3.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax3.set_ylabel('Temperature (K)',size=sfont)
ax3.set_title('Colors=HFE (MJ/kg)',size=sfont)
ax3.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax3)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig2.colorbar(im1, cax=cax, ticks = [0,.333,.6667,  1], orientation='vertical')
cbar.ax.set_yticklabels(['-1000','-500','0','500'])  # vertical colorbar
cbar.set_label('Log$_{10}$(HFE) [MJ/kg]',size=sfont)



#================ KPA flag
cmap2 = plt.get_cmap('plasma', 6)    # 6 discrete colors

im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.KPA)-0.5)/6.,cmap=cmap2,vmin=0,vmax=1,edgecolors='face')
#im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.KPA))/6.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax4.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax4.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax4.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax4.set_ylim(0.,20000)
ax4.set_xlim(0,30.)
ax4.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax4.set_ylabel('Temperature (K)',size=sfont)
ax4.set_title('Colors=KPA flag (1-6)',size=sfont)
ax4.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax4)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['0','1','2','3','4','5',  '6'])  # vertical colorbar
cbar.set_label('KPA Flag',size=sfont)


ax5.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.KPA)-0.5)/6.,cmap=cmap2,vmin=0,vmax=1)
ax5.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax5.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax5.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax5.set_ylim(10.,1.E7)
ax5.set_xlim(1.E-21,100.)
ax5.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax5.set_ylabel('Temperature (K)',size=sfont)
ax5.set_title('Colors=KPA flag (1-6)',size=sfont)
ax5.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax5)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
#cbar = fig.colorbar(im1, cax=cax, ticks = np.arange(13)/12, orientation='vertical')
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['','1','','2','','3','','4','','5','',  '6'])  # vertical colorbar
cbar.set_label('KPA Flag',size=sfont)
No description has been provided for this image
No description has been provided for this image

Model Data Quality (MDQ) Flag¶

We provide an MDQ rating which is saved as a 301-style SESAME table. The MDQ rating is a qualitative guide for users that provides some known quality control information. This information is necessarily subjective and incomplete. The rating does not guarantee any fidelity in the EOS model except for points that have been directly compared to experimental data in the plots below.

  • MDQ=1: Data point has experimental data and falls within experimental error. NOT USED AT THIS TIME. Will be added in the future.
  • MDQ=0: Data point rating is neutral. The value appears to be a thermodynamically self-consistent model point without known gross errors.
  • MDQ=-1: Data point rating is erroneous. The value is not thermodynamically correct in a known way. In this work, this flag is used for the replacement values for the tension region to generate an all positive pressure equilibrium EOS table. The current replacement scheme is to use the values in the adjacent grid cell: for each P and T in the sublimation region, use values at the adjacent density with positive pressure. The magnitude of the error is small, but the derivatives in this region will be incorrect.
  • MDQ=-2: Data point rating is non-physical. The value is not thermodynamically correct in a known way. In this work, this flag is used for ANEOS output with negative values for the specific entropy.
  • MDQ=-4: Developer judgement that this region of phase space does not provide a robust representation of the material. There are no high-pressure solid phases in this EOS model. Seek alternative models for problems focused in this region.
  • MDQ=-5: Developer judgement that this region of phase space does not provide a robust representation of the material. Applied to regions with pressures > 2 TPa. The Gruneisen gamma value is too large for liquid in this region and missing physics at higher pressures and temperatures. Seek alternative models for problems focused in this region.
In [13]:
# check the user-defined model data quality flag across the EOS surface
uniqmdqarr = np.unique(NewEOS.MDQ)
print("Unique values in the MDQ flag array:", uniqmdqarr)
print("Value: # elements with this value")
#get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x == y]

for i in uniqmdqarr:
    ix = get_indexes(i,np.reshape(NewEOS.MDQ,NewEOS.ND*NewEOS.NT))
    print(i,': ',len(ix))


##============= start new plot
fig3, ([ax4, ax5])= plt.subplots(nrows=1,ncols=2,figsize=(14,8))
plt.subplots_adjust(wspace=0.4)

#================ MDQ flag
cmap2 = plt.get_cmap('plasma', 6)    # 6 discrete colors

#im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=cmap2,vmin=0,vmax=1,edgecolors='face')
ax4.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax4.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax4.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax4.set_ylim(0.,20000)
ax4.set_xlim(0,30.)
ax4.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax4.set_ylabel('Temperature (K)',size=sfont)
ax4.set_title('Colors=MDQ flag',size=sfont)
ax4.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax4)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig3.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['','-5','','-4','','-3','','-2','','-1','','0','','1'])  # vertical colorbar

cmap2 = plt.get_cmap('plasma', 6)    # 6 discrete colors
ax5.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=cmap2,vmin=0,vmax=1)
#ax5.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax5.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax5.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax5.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax5.set_ylim(10.,1.E7)
ax5.set_xlim(1.E-21,100.)
ax5.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax5.set_ylabel('Temperature (K)',size=sfont)
ax5.set_title('Colors=MDQ flag',size=sfont)
ax5.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax5)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig3.colorbar(im1, cax=cax, ticks = np.arange(13)/12, orientation='vertical')
#cbar.ax.set_yticklabels(['','-5','','-4','','-3','','-2','','-1','','0','','1'])  # vertical colorbar
cbar.ax.set_yticklabels(['','-5','','-4','','-3','','-2','','-1','','0',''])  # vertical colorbar
Unique values in the MDQ flag array: [-2. -1.  0.]
Value: # elements with this value
-2.0 :  29077
-1.0 :  5726
0.0 :  530874
Out[13]:
[Text(1, 0.0, ''),
 Text(1, 0.08333333333333333, '-5'),
 Text(1, 0.16666666666666666, ''),
 Text(1, 0.25, '-4'),
 Text(1, 0.3333333333333333, ''),
 Text(1, 0.4166666666666667, '-3'),
 Text(1, 0.5, ''),
 Text(1, 0.5833333333333334, '-2'),
 Text(1, 0.6666666666666666, ''),
 Text(1, 0.75, '-1'),
 Text(1, 0.8333333333333334, ''),
 Text(1, 0.9166666666666666, '0'),
 Text(1, 1.0, '')]
No description has been provided for this image

ANEOS Model Compared to Experimental and Ab Initio Data¶

The following plots compare the ANEOS model Hugoniots with laboratory data and ab initio calculations.

ANEOS STP Hugoniots: blue -- Hugoniot based at STP (model solid at reference density).

In [14]:
# Display the initial conditions for the Hugoniot in ANEOS.OUTPUT    
print('\nRequested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):')
print('Initial Pressure (GPa):',NewEOS.anhug.P[0])
print('Initial Temperature (K):',NewEOS.anhug.T[0])
print('Initial Sp. Entropy (MJ/K/kg):',NewEOS.anhug.S[0])
print('Initial Sp. Internal Energy (MJ/kg):',NewEOS.anhug.U[0])
print('Initial Density (g/cm3):',NewEOS.anhug.rho[0])
#

#========== PLOTTING STYLES FOR LOTS OF PLOTS =================
styleqmd = ['^','red',5,'Fo QMD calc (Root et al. 2018)']
stylezdat = ['^','black',5,'Fo Z data (Root et al. 2018)']
stylegg = ['+','green',5,'Fo Gas gun data\n(Mosenfelder et al. 2007)']
styleliq = ['^','orange',7,'Liquid Fo data corrected\n(Thomas & Asimow 2013)']
styleliqhug = ['-','orange',1,'ANEOS Liquid 2273K Hugoniot']
stylezfit = ['-','green',1,'Fo Z fit (Davies et al. 2019)']
styleanhug = ['-','red',1,'ANEOS Glass-tension Hugoniot']
stylestphug = ['-','blue',1,'ANEOS Qtz Hugoniot']
styleanphase = ['-','black',1,'ANEOS Phase Boundaries']
stylephasept = ['o','black',6]



#======================LOTS OF PLOTS============================
# plot the Hugoniots for comparison
sfont = 15
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(14,18))
plt.subplots_adjust(wspace=0.4)

#-------------
# Pressure-Density Data vs. ANEOS
ai=0
aj=0
#axes[ai,aj].plot(NewEOS.anhug.rho,NewEOS.anhug.P,styleanhug[0],linewidth=styleanhug[2],color=styleanhug[1],label=styleanhug[3])

#axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
#axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    ind = np.where(qtz_hug_KD.P>50)[0]
    axes[ai,aj].plot(qtz_hug_KD.rho[ind],qtz_hug_KD.P[ind],'g-',linewidth=3,label='KD90 liq fit')
    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].rho,shockdat_list[L83qtzhug['id']].P,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].P_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].rho,shockdat_list[W62qtz['id']].P,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[A23qtz['id']].rho,shockdat_list[A23qtz['id']].P,fmt='gD',mfc='white',label=A23qtz['lab'],\
    #            xerr=shockdat_list[A23qtz['id']].rho_std,yerr=shockdat_list[A23qtz['id']].P_std)
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].rho,shockdat_list[T20qtz['id']].P,fmt='<',mfc='white',color='grey',\
                 xerr = shockdat_list[T20qtz['id']].rho_std, yerr=shockdat_list[T20qtz['id']].P_std, label=T20qtz['lab'])

    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].rho,shockdat_list[K13qtz['id']].P,fmt='gd',markersize=3,\
             xerr=shockdat_list[K13qtz['id']].rho_std,yerr=shockdat_list[K13qtz['id']].P_std,label=K13qtz['lab'])
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].rho,shockdat_list[L03qtz['id']].P,fmt='r>',label=L03qtz['lab'])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

axes[ai,aj].set_ylim(0,200)
axes[ai,aj].set_xlim(4.,5.5)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure - Temperature Data vs. ANEOS
ai=0
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T,'-.',color='orange',label='ANEOS FS Hugoniot')

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

if 0:
    axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T/1.e3,color='green',label='Qtz Brygoo et al. 2015')
    axes[ai,aj].plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,'-',color='red',label=fs_hug_Hicks06.label)
    axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T/1.E3,marker='+',color='green',label='Qtz Lyzenga & Ahrens 1980')
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T/1.E3,color='cyan',label=stv_mc_MM15.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[B88qtz['id']].P,shockdat_list[B88qtz['id']].T,fmt='go',mfc='white',yerr=shockdat_list[B88qtz['id']].T_std,label=B88qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].P,shockdat_list[L04qtz['id']].T,fmt='gv',mfc='white',xerr=shockdat_list[L04qtz['id']].P_std,yerr=shockdat_list[L04qtz['id']].T_std,label=L04qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[S81qtz['id']].P,shockdat_list[S81qtz['id']].T,fmt='gx',label=S81qtz['lab'])
    # Hicks decaying shock
    axes[ai,aj].plot(qtz_hug_Hicks06.P,qtz_hug_Hicks06.T,'g+:',label='H06 qtz')
    ind = np.where((shockdat_list[B88qtz['id']].P >75) & (shockdat_list[B88qtz['id']].P <110))
    qtztpfit = np.polyfit(shockdat_list[B88qtz['id']].P[ind],shockdat_list[B88qtz['id']].T[ind],2)
    # Hicks decaying shock
    #plt.plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T,'b+:',label='H06 FS')
    axes[ai,aj].errorbar(fs_hug_Hicks06.P,fs_hug_Hicks06.T,fmt='b+:',\
                 yerr=fs_hug_Hicks06.T_std,label='H06 FS')
    print(shockdat_list[B88qtz['id']].P[ind],shockdat_list[B88qtz['id']].T[ind])
    print('quartz Boslough P-T fit 75-110 GPa: ',qtztpfit)
    #plt.plot(shockdat_list[B88qtz['id']].P[ind],qtztpfit[1]+qtztpfit[0]*shockdat_list[B88qtz['id']].P[ind],'g:',label='Likely d-NiAs')
    qtzptmp = np.arange(75,110,1)
    qtzttmp = qtztpfit[2]+qtztpfit[1]*qtzptmp+qtztpfit[0]*qtzptmp*qtzptmp
    axes[ai,aj].plot(qtzptmp,qtzttmp/1.E3,'g--',label='CaCl2 region P-T fit')        
    axes[ai,aj].errorbar(shockdat_list[A23qtz['id']].P,shockdat_list[A23qtz['id']].T,fmt='gD',mfc='white',xerr=shockdat_list[A23qtz['id']].P_std,yerr=shockdat_list[A23qtz['id']].T_std,label=A23qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[B88fs['id']].P,shockdat_list[B88fs['id']].T,fmt='bo',mfc='white',yerr=shockdat_list[B88fs['id']].T_std,label=B88fs['lab'])

    axes[ai,aj].errorbar(shockdat_list[S80fs['id']].P,shockdat_list[S80fs['id']].T,fmt='bv',xerr=shockdat_list[S80fs['id']].P_std,yerr=shockdat_list[S80fs['id']].T_std,label=S80fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[S81fs['id']].P,shockdat_list[S81fs['id']].T,fmt='bx',label=S81fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].P,shockdat_list[L03qtz['id']].T,fmt='r>',label=L03qtz['lab'])

    
axes[ai,aj].set_xlim(0,200)
axes[ai,aj].set_ylim(0,10000)
axes[ai,aj].set_ylabel('Temperature (1000 K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Particle Velocity vs. Shock Velocity; Data vs. ANEOS
ai=1
aj=0

axes[ai,aj].plot(NewEOS.hug.up,NewEOS.hug.us,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.up,NewEOS.hugo.us,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.up,qtz_hug_W62.us,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.up,qtz_hug_L3.us,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.up,qtz_hug_KD.us,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.up,fs_hug_RTK19.us,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    mqtzup = deepcopy(np.asarray(qtzdf['up'][indqtz]))
    mqtzus = deepcopy(np.asarray(qtzdf['us'][indqtz]))
    #inda = np.where((mqtzus>8) & (mqtzus<9.23))[0]
    inda = np.where((mqtzus>6))[0]
    axes[ai,aj].errorbar(mqtzup[inda],mqtzus[inda],fmt='g+',label='M80 qtz included')

    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].up,shockdat_list[K13qtz['id']].us,fmt='gd',markersize=6,\
                 xerr=shockdat_list[K13qtz['id']].up_std,yerr=shockdat_list[K13qtz['id']].us_std,label=K13qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62fs['id']].up,shockdat_list[W62fs['id']].us,fmt='bx',label=W62fs['lab'])

    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].up,shockdat_list[L83qtzhug['id']].us,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].us_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].up,shockdat_list[W62qtz['id']].us,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[L03qtz['id']].rho,shockdat_list[L03qtz['id']].P,fmt='r>',label=L03qtz['lab'])
    #plt.errorbar(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,fmt='g<',mfc='white',\
    #             xerr = shockdat_list[T20qtz['id']].up_std, yerr=shockdat_list[T20qtz['id']].us_std, label=T20qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,fmt='<',mfc='white',color='grey',\
                  label=T20qtz['lab'])
    #print(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us)
    #res = np.polyfit(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,1)
    #print(res)
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].up,shockdat_list[L04qtz['id']].us,fmt='gv',mfc='white',\
             xerr=shockdat_list[L04qtz['id']].up_std,yerr=shockdat_list[L04qtz['id']].us_std,label=L04qtz['lab'])
    axes[ai,aj].errorbar(fs_hug_Hicks06.up,fs_hug_Hicks06.us,fmt='b+:',\
                label='H06 FS')
    axes[ai,aj].errorbar(shockdat_list[B19fs['id']].up,shockdat_list[B19fs['id']].us,fmt='b>',mfc='white',label=B19fs['lab'],\
                 xerr=shockdat_list[B19fs['id']].up_std,yerr=shockdat_list[B19fs['id']].us_std)
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].up,shockdat_list[L03qtz['id']].us,fmt='r>',label=L03qtz['lab'])

axes[ai,aj].set_xlim(0,6)
axes[ai,aj].set_ylim(0,12)
axes[ai,aj].set_xlabel('Particle Velocity (km/s)',size=sfont)
axes[ai,aj].set_ylabel('Shock Velocity (km/s)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Sp. Entropy vs Pressure; Calculated from experiments vs. ANEOS
ai=1
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.S*1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.S*1.e3,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    # Kraus entropy fit only good to 800 GPa
    ind = np.where(qtz_hug_KD.P <= 800.)
    axes[ai,aj].plot(qtz_hug_KD.P[ind],qtz_hug_KD.S[ind],color='green',label='Qtz Kraus et al. 2012')
    Kraus_serr = (127.3+0.328*qtz_hug_KD.P[ind]-0.000193*np.power(qtz_hug_KD.P[ind],2))/1.e3 # kJ/K/kg
    axes[ai,aj].fill_between(qtz_hug_KD.P[ind],(qtz_hug_KD.S[ind]-Kraus_serr),(qtz_hug_KD.S[ind]+Kraus_serr),color=stylezfit[1],alpha=0.2)

# silicadb
if 1:
    entropyshift = 220.
    penter = 111.
    kqtzparr = np.arange(100,500) # GPa
    kqtzsarr = 1793+22.4*kqtzparr-0.0426*np.power(kqtzparr,2)+4.319e-05*np.power(kqtzparr,3)-1.742e-8*np.power(kqtzparr,4)
    kqtzserrarr = 127.3+0.328*kqtzparr-0.000193*np.power(kqtzparr,2)
    ind = np.where(kqtzparr >= penter)[0]
    axes[ai,aj].plot(kqtzparr[ind],kqtzsarr[ind]/1.e3,'g-',linewidth=1,label='K12 qtz Hugoniot-liq region')
    axes[ai,aj].plot(kqtzparr[ind],(kqtzsarr[ind]-entropyshift)/1.e3,'g-',linewidth=3,label='K12 Qtz Hugoniot shift')
    axes[ai,aj].fill_between(kqtzparr[ind],(kqtzsarr[ind]-entropyshift-kqtzserrarr[ind])/1.e3,(kqtzsarr[ind]-entropyshift+kqtzserrarr[ind])/1.e3,color=stylezfit[1],alpha=0.2)
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].P,shockdat_list[L03qtz['id']].S,fmt='r>',label=L03qtz['lab'])
   
    

axes[ai,aj].set_ylim(0,8.)
axes[ai,aj].set_xlim(0,250)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   LOW PRESSURE + LIQUID HUGONIOT
ai=2
aj=0
# ANEOS REFERENCE POINT (STP) Hugoniots
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].rho,shockdat_list[L83qtzhug['id']].P,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].P_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].rho,shockdat_list[W62qtz['id']].P,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[A23qtz['id']].rho,shockdat_list[A23qtz['id']].P,fmt='gD',mfc='white',label=A23qtz['lab'],\
    #            xerr=shockdat_list[A23qtz['id']].rho_std,yerr=shockdat_list[A23qtz['id']].P_std)
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].rho,shockdat_list[T20qtz['id']].P,fmt='<',mfc='white',color='grey',\
                 xerr = shockdat_list[T20qtz['id']].rho_std, yerr=shockdat_list[T20qtz['id']].P_std, label=T20qtz['lab'])
    #print(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us)
    #res = np.polyfit(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,1)
    #print(res)
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].rho,shockdat_list[L04qtz['id']].P,fmt='gv',mfc='white',\
             xerr=shockdat_list[L04qtz['id']].rho_std,yerr=shockdat_list[L04qtz['id']].P_std,label=L04qtz['lab'])

    ind = np.where(qtz_hug_KD.P>112)[0]
    axes[ai,aj].plot(qtz_hug_KD.rho[ind],qtz_hug_KD.P[ind],'g-',linewidth=3,label='KD90 liq fit')
    axes[ai,aj].errorbar(qtzdf['rho'][indqtz],qtzdf['P'][indqtz],fmt='+',color='grey',label='M80 qtz')

    
axes[ai,aj].set_ylim(-5,200)
axes[ai,aj].set_xlim(2.5,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   EXTRAPOLATED HIGH PRESSURE
ai=2
aj=1
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[W62fs['id']].rho,shockdat_list[W62fs['id']].P,fmt='bx',label=W62fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[L83fshug['id']].rho,shockdat_list[L83fshug['id']].P,fmt='bo',label=L83fshug['lab'])
    axes[ai,aj].errorbar(shockdat_list[B19fs['id']].rho,shockdat_list[B19fs['id']].P,fmt='b>',xerr=shockdat_list[B19fs['id']].rho_std,\
                 yerr=shockdat_list[B19fs['id']].P_std,mfc='white',label=B19fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[T18fs['id']].rho,shockdat_list[T18fs['id']].P,fmt='b<',mfc='white',\
                 xerr=shockdat_list[T18fs['id']].rho_std,yerr=shockdat_list[T18fs['id']].P_std,label=T18fs['lab'])

    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho)/1000.,shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'])
    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho),shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'],\
    #            xerr=shockdat_list[A23fs['id']].rho_std,yerr=shockdat_list[A23fs['id']].P_std)
    axes[ai,aj].errorbar(fsdf['rho'][indfs],fsdf['P'][indfs],fmt='+',color='grey',label='M80 FS excluded')
    indhigh = np.where(fsdf['P'][indfs] > 74)[0]
    axes[ai,aj].errorbar(fsdf['rho'][indfs][indhigh],fsdf['P'][indfs][indhigh],fmt='b+',label='M80 FS included')

    #axes[ai,aj].plot(fitfslrarr,fitfslparr,'b-',linewidth=2,label='cubic fit liquid')
    #axes[ai,aj].plot(fitfssrarr,fitfssparr,'b-',label='linear fit solid')
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)


axes[ai,aj].set_ylim(0,300)
axes[ai,aj].set_xlim(3.,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

plt.savefig(MODELNAME+'-Hugoniots.pdf')
Requested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):
Initial Pressure (GPa): 0.0001
Initial Temperature (K): 298.0
Initial Sp. Entropy (MJ/K/kg): 0.00017886
Initial Sp. Internal Energy (MJ/kg): 0.043369
Initial Density (g/cm3): 2.2
[ 75.9  85.9  92.5  99.3 107.8 109.7] [4440. 4850. 5430. 5630. 5900. 6100.]
quartz Boslough P-T fit 75-110 GPa:  [-3.39692149e-01  1.11592149e+02 -2.10423906e+03]
No description has been provided for this image
In [15]:
print(shockdat_list[L03qtz['id']].S)
[3.365 3.527 3.621 3.712 3.817 3.839]
In [16]:
# Display the initial conditions for the Hugoniot in ANEOS.OUTPUT    
print('\nRequested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):')
print('Initial Pressure (GPa):',NewEOS.anhug.P[0])
print('Initial Temperature (K):',NewEOS.anhug.T[0])
print('Initial Sp. Entropy (MJ/K/kg):',NewEOS.anhug.S[0])
print('Initial Sp. Internal Energy (MJ/kg):',NewEOS.anhug.U[0])
print('Initial Density (g/cm3):',NewEOS.anhug.rho[0])
#

#========== PLOTTING STYLES FOR LOTS OF PLOTS =================
styleqmd = ['^','red',5,'Fo QMD calc (Root et al. 2018)']
stylezdat = ['^','black',5,'Fo Z data (Root et al. 2018)']
stylegg = ['+','green',5,'Fo Gas gun data\n(Mosenfelder et al. 2007)']
styleliq = ['^','orange',7,'Liquid Fo data corrected\n(Thomas & Asimow 2013)']
styleliqhug = ['-','orange',1,'ANEOS Liquid 2273K Hugoniot']
stylezfit = ['-','green',1,'Fo Z fit (Davies et al. 2019)']
styleanhug = ['-','red',1,'ANEOS Glass-tension Hugoniot']
stylestphug = ['-','blue',1,'Tabulated STP Hugoniot']
styleanphase = ['-','black',1,'ANEOS Phase Boundaries']
stylephasept = ['o','black',6]



#======================LOTS OF PLOTS============================
# plot the Hugoniots for comparison
sfont = 15
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16,12))
plt.subplots_adjust(wspace=0.4)

#-------------
# Pressure-Density Data vs. ANEOS
ai=0
aj=0
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)


axes[ai,aj].set_ylim(0,1400)
axes[ai,aj].set_xlim(3.,9.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure - Temperature Data vs. ANEOS
ai=0
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T/1.E3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label='ANEOS STP Hugoniot')
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T/1.E3,'-.',color='lime',label='ANEOS Glass Hugoniot')

if 0:
    axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T/1.e3,color='green',label='Qtz Brygoo et al. 2015')
    axes[ai,aj].plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,'-',color='red',label=fs_hug_Hicks06.label)

    axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T/1.E3,marker='+',color='green',label='Qtz Lyzenga & Ahrens 1980')
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T/1.E3,color='cyan',label=stv_mc_MM15.label)


axes[ai,aj].set_xlim(0,1200)
axes[ai,aj].set_ylim(0,70)
axes[ai,aj].set_ylabel('Temperature (1000 K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)


#-------
# pressure - temperature melt curve
ai=1
aj=0
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.P[0],NewEOS.hug.T[0],'o',color=stylestphug[1])

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,15000)
axes[ai,aj].set_xlim(-5,600)

#axes[ai,aj].plot(peridotiteliquidusdatarr.iloc[:,0].values,peridotiteliquidusdatarr.iloc[:,1].values,linestyle='-',color='orange',label='Peridotite Liquidus')
#axes[ai,aj].plot(peridotitesolidusdatarr.iloc[:,0].values,peridotitesolidusdatarr.iloc[:,1].values,linestyle=':',color='orange',label='Peridotite Solidus\n(Fiquet et al. 2010)')

if 1:
    #axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T,color='cyan',label=stv_mc_MM15.label)
    axes[ai,aj].plot(qtz_mc_BR98.Pl,qtz_mc_BR98.T,color='pink',label=qtz_mc_BR98.label)

if 1:
    wcol = "#653700"
    wcol2 = "#FBDD7E"
    wcol3 = "#F97306"
    wop = 0.5
    wop2=0.75
    wop3=0.25
    ind = np.arange(0,11)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],label='melt curve exp',alpha=wop,color=wcol)
    ind = np.arange(12,16)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
    ind = np.arange(17,21)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
    ind = np.arange(22,24)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2,label='melt curve calc')
    ind = np.arange(25,28)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2)
    # Millot stv melt curve crossing
    axes[ai,aj].fill_between([470.,530.],[8000.,8000.],[8600.,8600.],alpha=0.5,color='blue',label='MM15')
    

if 0:
    ind = np.where(stsmcsmooth.P <=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m-',linewidth=3,label='Proposed melt curve')
    ind = np.where(stsmcsmooth.P >=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m:',linewidth=3)#,label='Proposed melt curve')



axes[ai,aj].legend()

#-------
ai=1
aj=1
#axes[ai,aj].axis('off')

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.S*1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.S*1.e3,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    # Kraus entropy fit only good to 800 GPa
    ind = np.where(qtz_hug_KD.P <= 800.)
    axes[ai,aj].plot(qtz_hug_KD.P[ind],qtz_hug_KD.S[ind],color='green',label='Qtz Kraus et al. 2012')
    Kraus_serr = (127.3+0.328*qtz_hug_KD.P[ind]-0.000193*np.power(qtz_hug_KD.P[ind],2))/1.e3 # kJ/K/kg
    axes[ai,aj].fill_between(qtz_hug_KD.P[ind],(qtz_hug_KD.S[ind]-Kraus_serr),(qtz_hug_KD.S[ind]+Kraus_serr),color=stylezfit[1],alpha=0.2)

# silicadb
if 1:
    entropyshift = 220.
    penter = 111.
    kqtzparr = np.arange(100,500) # GPa
    kqtzsarr = 1793+22.4*kqtzparr-0.0426*np.power(kqtzparr,2)+4.319e-05*np.power(kqtzparr,3)-1.742e-8*np.power(kqtzparr,4)
    kqtzserrarr = 127.3+0.328*kqtzparr-0.000193*np.power(kqtzparr,2)
    ind = np.where(kqtzparr >= penter)[0]
    axes[ai,aj].plot(kqtzparr[ind],kqtzsarr[ind]/1.e3,'g-',linewidth=1,label='K12 qtz Hugoniot-liq region')
    axes[ai,aj].plot(kqtzparr[ind],(kqtzsarr[ind]-entropyshift)/1.e3,'g-',linewidth=3,label='K12 Qtz Hugoniot shift')
    axes[ai,aj].fill_between(kqtzparr[ind],(kqtzsarr[ind]-entropyshift-kqtzserrarr[ind])/1.e3,(kqtzsarr[ind]-entropyshift+kqtzserrarr[ind])/1.e3,color=stylezfit[1],alpha=0.2)
   
    

axes[ai,aj].set_ylim(0,5.)
axes[ai,aj].set_xlim(0,200)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

plt.savefig(MODELNAME+'-3panel.pdf')
Requested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):
Initial Pressure (GPa): 0.0001
Initial Temperature (K): 298.0
Initial Sp. Entropy (MJ/K/kg): 0.00017886
Initial Sp. Internal Energy (MJ/kg): 0.043369
Initial Density (g/cm3): 2.2
[13.7  14.5  20.   23.75 28.   34.   36.   37.5  40.   42.   46.  ]
[55. 60. 66. 74.]
[100. 102. 114. 117.]
[135. 155.]
[217.5 218.5   nan]
No description has been provided for this image

Model Critical Point and Shock-Induced Phase Changes¶

In [17]:
print('ANEOS Calculated Critical point:')
print('   Temperature (K):       ',NewEOS.cp.T)
print('   Pressure (GPa):        ', NewEOS.cp.P)
print('   Density (g/cm^3):      ',NewEOS.cp.rho)
print('   Sp. Energy (MJ/kg):    ', NewEOS.cp.U)
print('   Sp. Entropy (kJ/K/kg): ', NewEOS.cp.S*1.e3)
print('\n')
print('ANEOS Calculated Values at 1 bar:')
print('   Temperature at melting (K):                     ', NewEOS.onebar.Tmelt)
print('   Liquid density at melting (g/cm^3):             ', NewEOS.onebar.rhocm)
print('   Sp. Entropy at inc. melting Sim (kJ/K/kg):      ', NewEOS.onebar.Sim*1.E3)
print('   Sp. Entropy at com. melting Scm (kJ/K/kg):      ', NewEOS.onebar.Scm*1.E3)
print('   Temperature at boiling (K):                     ', NewEOS.onebar.Tvap)
print('   Liquid density at boiling (g/cm^3):             ', NewEOS.onebar.rhoiv)
print('   Sp. Entropy at inc. vaporization Siv (kJ/K/kg): ', NewEOS.onebar.Siv*1.E3)
print('   Sp. Entropy at com. vaporization Scv (kJ/K/kg): ', NewEOS.onebar.Scv*1.E3)
print('\n')
if NewEOS.tp.T > 0:
    print('ANEOS Calculated Triple Point:')
    print('   Temperature (K):                        ',NewEOS.tp.T)
    print('   Pressure (GPa):                         ', NewEOS.tp.P)
    print('   Solid density (g/cm^3):                 ',NewEOS.tp.rhos)
    print('   Liquid density (g/cm^3):                ',NewEOS.tp.rhol)
    print('   Vapor density (g/cm^3):                 ',NewEOS.tp.rhov)
    print('   Sp. Entropy at inc. melt Sim (kJ/K/kg): ', NewEOS.tp.Sim*1.E3)
    print('   Sp. Entropy at com. melt Scm (kJ/K/kg): ', NewEOS.tp.Scm*1.E3)
    print('   Sp. Entropy at inc. vap. Siv (kJ/K/kg): ', NewEOS.tp.Siv*1.E3)
    print('   Sp. Entropy at com. vap. Scv (kJ/K/kg): ', NewEOS.tp.Scv*1.E3)
else:
    print('ANEOS: No melt curve or triple point in this calculation.')
print('\n')
print('Starting at reference conditions (rho,T): ',NewEOS.R0REF,NewEOS.T0REF, ' (g/cm^3, K)')
print('ANEOS model conditions for shock-induced vaporization on release to 1 bar:')
print('                     :         IM        CM        IV      50%V        CV')
print('Shock Pressure (GPa) : {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.onebar.Sim,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Scm,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Siv,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Siv+0.5*(NewEOS.onebar.Scv-NewEOS.onebar.Siv),NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Scv,NewEOS.hug.S,NewEOS.hug.P)))
print('Imp. Match Vel (km/s): {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.onebar.Sim,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Scm,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Siv,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Siv+0.5*(NewEOS.onebar.Scv-NewEOS.onebar.Siv),NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Scv,NewEOS.hug.S,NewEOS.hug.up)*2))
print('ANEOS model conditions for shock-induced vaporization on release to the triple point pressure:')
print('                     :         IM        CM        IV      50%V        CV')
print('Shock Pressure (GPa) : {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.tp.Sim,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Scm,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Siv,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Siv+0.5*(NewEOS.tp.Scv-NewEOS.tp.Siv),NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Scv,NewEOS.hug.S,NewEOS.hug.P)))
print('Imp. Match Vel (km/s): {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.tp.Sim,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Scm,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Siv,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Siv+0.5*(NewEOS.tp.Scv-NewEOS.tp.Siv),NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Scv,NewEOS.hug.S,NewEOS.hug.up)*2))
ANEOS Calculated Critical point:
   Temperature (K):        7512.49
   Pressure (GPa):         0.390284
   Density (g/cm^3):       0.507337
   Sp. Energy (MJ/kg):     16.8021
   Sp. Entropy (kJ/K/kg):  5.76996


ANEOS Calculated Values at 1 bar:
   Temperature at melting (K):                      1996.010972480397
   Liquid density at melting (g/cm^3):              1.8810072665433093
   Sp. Entropy at inc. melting Sim (kJ/K/kg):       2.4828055952383483
   Sp. Entropy at com. melting Scm (kJ/K/kg):       2.8135042388169302
   Temperature at boiling (K):                      3212.86
   Liquid density at boiling (g/cm^3):              1.6832
   Sp. Entropy at inc. vaporization Siv (kJ/K/kg):  3.58174
   Sp. Entropy at com. vaporization Scv (kJ/K/kg):  6.80183


ANEOS Calculated Triple Point:
   Temperature (K):                         1996.0
   Pressure (GPa):                          2.90104e-09
   Solid density (g/cm^3):                  2.0225999999999997
   Liquid density (g/cm^3):                 1.881
   Vapor density (g/cm^3):                  7.5663e-09
   Sp. Entropy at inc. melt Sim (kJ/K/kg):  2.4828
   Sp. Entropy at com. melt Scm (kJ/K/kg):  2.8135
   Sp. Entropy at inc. vap. Siv (kJ/K/kg):  2.81349
   Sp. Entropy at com. vap. Scv (kJ/K/kg):  8.45833


Starting at reference conditions (rho,T):  2.2 298.0  (g/cm^3, K)
ANEOS model conditions for shock-induced vaporization on release to 1 bar:
                     :         IM        CM        IV      50%V        CV
Shock Pressure (GPa) :      39.70     50.46     84.62    243.49    725.55
Imp. Match Vel (km/s):       6.05      6.86      9.01     15.92     28.78
ANEOS model conditions for shock-induced vaporization on release to the triple point pressure:
                     :         IM        CM        IV      50%V        CV
Shock Pressure (GPa) :      39.70     50.46     50.46    330.10   1872.74
Imp. Match Vel (km/s):       6.05      6.86      6.86     18.79     47.50

Model Phase Boundaries¶

Black lines are the ANEOS phase boundaries. Blue curve is the interpolated Hugoniot from the gridded eos table.

Comparison critical points from experimentally constrained quartz critical point (Kraus et al. 2012), ab initio enstatite critical point (Xiao & Stixrude 2018), ab initio forsterite critical point (Townsend et al. 2020).

Orange diamonds are the 1-bar boiling point; orange squares are the 1-bar melting point for pure forsterite for reference.

In [18]:
# PHASE DIAGRAM PLOTS
#======================LOTS OF PLOTS============================
# Plot aneos phase boundaries with focus on the vapor curve
sfont = 15
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(14,24))
plt.subplots_adjust(wspace=0.4)

#------------------------------
# entropy-temperature
ai=0
aj=0
axes[ai,aj].plot(NewEOS.vc.Sl*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Sv*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.S*1.E3,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Ss*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.S*1.E3,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.S[0]*1.E3,NewEOS.hug.T[0],'o',color=stylestphug[1])

axes[ai,aj].set_ylim(0,1.5*NewEOS.cp.T)
axes[ai,aj].set_xlim(0,14)
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Temperature [K]',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------------------
# density-pressure
ai=0
aj=1
axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.rho[0],NewEOS.hug.P[0],'o',color=stylestphug[1])



axes[ai,aj].set_ylim(1.e-5,300)
axes[ai,aj].set_xlim(1.e-1,8)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#------
#entropy-pressure
ai=1
aj=0
axes[ai,aj].semilogy(NewEOS.vc.Sl*1.E3,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Sv*1.E3,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.S*1.E3,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.E3,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Ss*1.E3,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.S*1.E3,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.S[0]*1.E3,NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(1.e-10,10000)
axes[ai,aj].set_xlim(0,14)

#-----------
# specific energy - pressure
ai=1
aj=1
axes[ai,aj].semilogy(NewEOS.vc.Ul,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Uv,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.U,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.U,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.U[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(1.e-11,10000)
axes[ai,aj].set_xlim(-3,NewEOS.cp.U*2.)

#-------
# temperature - pressure
ai=2
aj=0
axes[ai,aj].semilogy(NewEOS.vc.T,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.T,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.T,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.T,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.T,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.T,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.T[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Temperature (K)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_xlim(0,9000)
axes[ai,aj].set_ylim(1.e-16,100)

#-------
# pressure - temperature melt curve
ai=2
aj=1
#axes[ai,aj].plot(NewEOS.vc.Pl,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
#axes[ai,aj].plot(NewEOS.vc.Pv,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
#axes[ai,aj].plot(NewEOS.cp.P,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.P[0],NewEOS.hug.T[0],'o',color=stylestphug[1])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T,'-.',color='orange',label='ANEOS FS Hugoniot')

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,9000)
axes[ai,aj].set_xlim(-5,600)


#-------
# density - pressure melt curve
ai=3
aj=0
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.rho[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,300)
axes[ai,aj].set_xlim(2.5,6.5)


#-------
# density -- temperature
ai=3
aj=1
axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
#axes[ai,aj].plot(NewEOS.hug.T[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,9000)
axes[ai,aj].set_xlim(0.1,10)



    
#===================================================================
# PLOT EXPERIMENTAL CONSTRAINTS ON THE ANEOS PHASE BOUNDARIES
PLOTDATA=1
datms = 8
if PLOTDATA:
    # entropy-temperature [0,0]
    ai=0
    aj=0
    
    # Connolly JGR 2016 summary of critical points
    axes[ai,aj].add_patch(Rectangle((4.8,5071), 0.5, 1000,  alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')
    
    axes[ai,aj].plot(5.149,5130,'^',markersize=datms,color='black',label='Qtz Kraus et al. 2012')
    # Triple point from Nagahara et al. 1994
    #Triple point of forsterite 5.2e-5 bar, 1890 C (2163 K)
    axes[ai,aj].plot(3.470,2163,'X',color='orange',label='Fo TP Nagahara et al. 1994')

    axes[ai,aj].add_patch(Rectangle((1850, 1.e-10), 600, 3.e-7, alpha=0.2, facecolor='brown', linewidth=0))
    axes[ai,aj].plot([3.507,7.267],[3127,3127],'d',color='red',label='JANAF bp')
    axes[ai,aj].plot([3.507],[1996],'s',color='red',label='MK88 tp')

    
    if 1:
        #plot Connolly calcs
        b = np.argsort(C16highs)
        axes[ai,aj].plot(C16highs[b]/1.e3,C16hight[b],'-',color='black',label='C16 HSC-vWRK vapor curve')
        axes[ai,aj].plot(5.382,6074,'o',color='black',markersize=5)
        axes[ai,aj].plot(C16lows/1.e3,C16lowt,'--',color='black',label='C16 HSC-CSRK vapor curve')
        axes[ai,aj].plot(5.322,5548,'o',color='black',markersize=5)
        #plt.plot(C16sils,C16silt,'-',color='orangered',label='C16 SiO(l)-vWRK')
        #plt.plot(5403,5675,'o',color='orangered')
        # plot boiling point at 1 bar
        # note Kraus S_vap too small
        axes[ai,aj].plot([3.554,7.460],[3130,3130],'kD',label='C16 1-bar boiling pt.')
      
    # density-pressure [0,1]
    ai=0
    aj=1
    # Connolly JGR 2016 summary of critical points
    axes[ai,aj].add_patch(Rectangle((0.5,.13), 0.5, .8,  alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')

    axes[ai,aj].plot(0.4822,0.1414,'+',label='En Xiao & Stixrude 2018',markersize=datms,color='brown')
    axes[ai,aj].plot(0.508,0.132,'^',label='Qtz Kraus et al. 2012',markersize=datms,color='black')
    axes[ai,aj].plot(0.51,0.114,'*',label='Fo Townsend et al. 2020',markersize=datms,color='magenta')

    axes[ai,aj].plot(60.083/27.30,1.e-4,'s',color='red',label='Bourova Richet 1998 qtz')
    axes[ai,aj].plot(60.083/23.45,1.e-4,'d',color='orange',label='Bourova Richet 1998 liq')
    axes[ai,aj].plot(60.083/23.51,1.e-4,'d',color='orange',label='')

    # Caracas Pyrolite DFT
    #axes[ai,aj].add_patch(Rectangle((.3, .08), 0.7, 0.05, alpha=0.5, facecolor='cyan', linewidth=0))
    
    
    # entropy -- pressure [1,0]
    ai=1
    aj=0
    # 0.6689 kJ/kg STP entropy from Robie et al. 1982
    axes[ai,aj].plot(0.669,1.e-4,'+',color='orange',label='STP')
    # 1 bar melting point entropy, enthalpy and temperature from Richet et al. 1993
    axes[ai,aj].plot(3.470,1.e-4,'d',color='orange',label='MP 1 bar Richet et al. 1993')
    # Triple point from Nagahara et al. 1994
    #Triple point of forsterite 5.2e-5 bar, 1890 C (2163 K)
    axes[ai,aj].plot(3.470,5.2e-9,'X',color='orange',label='TP Nagahara et al. 1994')

    axes[ai,aj].plot([3.507,7.267],[1.e-4,1.e-4],'d',color='red',label='JANAF bp')

    
    # temperature -- pressure [2,0]
    ai=2
    aj=0
    
    # Petaev BSE triple system region
    #axes[ai,aj].add_patch(Rectangle((1850, 1.e-10), 600, 3.e-7, alpha=0.2, facecolor='brown', linewidth=0))
    #axes[ai,aj].scatter(1900, -1.e-9, alpha=0.2, marker='s',color='brown', label='Petaev BSE TP')
    
    # Triple point from Nagahara et al. 1994
    #Triple point of forsterite 5.2e-5 bar, 1890 C (2163 K)
    axes[ai,aj].plot(2163.,5.2e-9,'X',color='orange',label='Fo TP Nagahara et al. 1994')
    # Melt curve points
    # axes[ai,aj].plot(2174,1.e-4,'d',color='orange',label='Richet et al. 1993')
    # axes[ai,aj].plot(2550,13.,'v',color='orange',label='de Koker et al. 2009')
    # Nagahara et al. 1994 vapor curve of forsterite
    tvapfo = np.arange(46)*10.+1400.+273. # in C to K
    #pvapfo = np.exp(-77955.0/tvapfo+27.7)
    pvapfonagaqhara = np.exp(169.-543.e3/(8.314*(tvapfo))) 
    #print(pvapfonagahara)
    #print(tvapfo)
    #print(pvapfo)
    #axes[3,0].plot(tvapfo+273.,pvapfo*1.e5/1.e9,linestyle='-',color='orange',label='Sublimation curve\nMysen & Kushiro 1988')
    #axes[3,0].plot(tvapfo,pvapfonagahara*1.e5/1.e9,linestyle='-',color='green',label='Sublimation curve\nNagahara et al. 1994')
    #print(pvapfonagahara)
    #tt=np.zeros(7)
    #tt.fill(273.)
    #tt=tt+np.asarray([1400.,1500.,1600.,1700.,1725.,1760.,1800.])
    #axes[ai,aj].plot(tt,np.asarray([8.79E-9,5.98E-8,3.85E-7,2.7E-6,4.21E-6,9.58E-6,1.52E-5])/1.e4,'-x',color='orange',label='Fo Nagahara et al. 1994')

    #Centolanzi and Chapman Tektite vapor curve
    tvaptek = np.arange(1400)+1700. # in K
    pvaptek = np.exp(-57400.0/tvaptek+18.5)*101325. # in atm to Pa
    axes[ai,aj].plot(tvaptek,pvaptek/1.e9,linestyle='-',color='cyan',label='Vapor curve tektites\nCentolanzi & Chapman 1966')

    # Misha Petaev BSE vapor curve 50% condensed
    #axes[ai,aj].scatter([2500., 3000., 3500., 4000., 4500.],[0.5e-6,2.e-5,3.e-4,20e-4,100.e-4],marker='+',label='Petaev BSE 50% condensation')
    #axes[ai,aj].plot(6110,0.114,'*',label='Fo Townsend et al. 2020',markersize=datms,color='magenta')

    # Mysen and Kushiro 1988 SiO2 vapor curve and TPs
    # not sure this is pure SiO2 - may be mixed with H2 - check
    #axes[ai,aj].add_patch(Rectangle((1608+273, 4.e-11), 93, 1.3e-10, alpha=0.7, facecolor='green', linewidth=0))
    #axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='green', label='Mysen1988 TP')
    # fitted temps between 2000 and 1880 - this is SiO2+H2
    #MKvap_t = np.arange(1880,2000,20) # K Fig 5
    #MKvap_p = np.exp(-50011/MKvap_t + 14.2)*1.e5/1.e9 # GPa
    #axes[ai,aj].plot(MKvap_t,MKvap_p,color='red',label='Mysen1988 SiO2+H2 VC')
    #print(MKvap_t,MKvap_p)
    # Figure 2 pure SiO2 Hidalgo 1960 equation 1660-1920 K
    MKvap_t = np.arange(1660,1920,20) # K Fig 2
    MKvap_p = np.exp(-57800/MKvap_t + 18.4)*1.e5/1.e9 # GPa
    axes[ai,aj].plot(MKvap_t,MKvap_p,color='red',label='MK88/H60 SiO2 VC')
    
    # Connolly JGR 2016 summary of critical points
    axes[ai,aj].add_patch(Rectangle((5071, .13), 1000, .8, alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')
    
    # pressure- temperature [2,1]
    ai=2
    aj=1
    # Melt curve points
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T,label=stv_mc_MM15.label)
    axes[ai,aj].plot(qtz_mc_BR98.Pl,qtz_mc_BR98.T,label=qtz_mc_BR98.label)

    #axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T,label=qtz_hug_LA80.label)
    #axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T,color='green',label='Qtz Brygoo et al. 2015')


    #axes[ai,aj].plot(Mosenfelderdatarr.iloc[:,0].values,Mosenfelderdatarr.iloc[:,1].values,linestyle='-',color='brown',linewidth=1,label='Forsterite Melt Curve\n(Mosenfelder et al. 2007)')
    #axes[ai,aj].plot(peridotiteliquidusdatarr.iloc[:,0].values,peridotiteliquidusdatarr.iloc[:,1].values,linestyle='-',color='orange',label='Peridotite Liquidus')
    #axes[ai,aj].plot(peridotitesolidusdatarr.iloc[:,0].values,peridotitesolidusdatarr.iloc[:,1].values,linestyle=':',color='orange',label='Peridotite Solidus\n(Fiquet et al. 2010)')
    if 1:
        wcol = "#653700"
        wcol2 = "#FBDD7E"
        wcol3 = "#F97306"
        wop = 0.5
        wop2=0.75
        wop3=0.25
        ind = np.arange(0,11)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],label='melt curve exp',alpha=wop,color=wcol)
        ind = np.arange(12,16)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
        ind = np.arange(17,21)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
        ind = np.arange(22,24)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2,label='melt curve calc')
        ind = np.arange(25,28)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2)
        # Millot stv melt curve crossing
        axes[ai,aj].fill_between([470.,530.],[8000.,8000.],[8600.,8600.],alpha=0.5,color='blue',label='MM15')
 

    if 0:
        ind = np.where(stsmcsmooth.P <=Pt3)[0]
        axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m-',linewidth=3,label='Proposed melt curve')
        ind = np.where(stsmcsmooth.P >=Pt3)[0]
        axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m:',linewidth=3)#,label='Proposed melt curve')

    # density - pressure [3,0]
    ai=3
    aj=0
    
    #axes[ai,aj].scatter(qtz_hug_LA80.rho[6],qtz_hug_LA80.P[6],color='green',label=qtz_hug_LA80.label)
    #axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)

    # density - temperature [3,1]
    ai=3
    aj=1
    # Connolly JGR 2016 summary of critical points for silica
    axes[ai,aj].add_patch(Rectangle((0.5,5071), 0.5, 1000,  alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')

    axes[ai,aj].plot(60.083/27.30,1999.,'s',color='red',label='Bourova Richet 1998')
    axes[ai,aj].plot(60.083/23.45,1700.,'d',color='orange',label='Bourova Richet 1998')
    axes[ai,aj].plot(60.083/23.51,1673.,'d',color='orange',label='')

    #axes[ai,aj].scatter(qtz_hug_LA80.rho[6],qtz_hug_LA80.T[6],color='green',label=qtz_hug_LA80.label)

    # density - pressure [3,0]
    ai=3
    aj=1    
    # Connolly JGR 2016
    axes[ai,aj].plot(C16bsrho,C16bstemp,'k-',label='C16 HSC-vWRK')
    axes[ai,aj].plot(C16bdrho,C16bdtemp,'k--',label='C16 HSC-CSRK')
    axes[ai,aj].plot(C16rsrho,C16rstemp,'r-',label='C16 SiO(l)-vWRK')
    # critical points
    axes[ai,aj].plot(1.036,6074,'ko')
    axes[ai,aj].plot(.555,5548,'ko')
    axes[ai,aj].plot(.782,5675,'ro')

axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)
axes[2,0].legend(fontsize=sfont-5)
axes[2,1].legend(fontsize=sfont-5)
axes[3,0].legend(fontsize=sfont-5)
axes[3,1].legend(fontsize=sfont-5)




plt.tight_layout(pad=0.6)
plt.savefig(MODELNAME+'-phase-diagram-fig.pdf', format='pdf', dpi=1000,transparent=True)
#--------
[13.7  14.5  20.   23.75 28.   34.   36.   37.5  40.   42.   46.  ]
[55. 60. 66. 74.]
[100. 102. 114. 117.]
[135. 155.]
[217.5 218.5   nan]
No description has been provided for this image

ANEOS model comparisons to experimental data¶

In the liquid field, we assume that the heat capacity is constant and similar to the value for forsterite liquid reported in Thomas and Asimow 2013.

In [19]:
# 1-bar plots and 298 K plots

# Plot aneos vapor curves
sfont = 15
fig, axes = plt.subplots(2, 2, figsize=(14,14))
plt.subplots_adjust(wspace=0.4)

#--------
# temperature - entropy at 1 bar
ai=0
aj=0

axes[ai,aj].plot(NewEOS.onebar.T,NewEOS.onebar.S*1.E3,'+',color='blue',label='ANEOS 1 bar')

axes[ai,aj].plot(qtz_1bar_JANAF.T,qtz_1bar_JANAF.S,linestyle='-',color='brown',label=qtz_1bar_JANAF.label)
#axes[ai,aj].plot(janafgas_1bar.T,janafgas_1bar.S/1.e3,linestyle='-',color='orange',label='JANAF gas')
axes[ai,aj].scatter([3127,3127],[3.507,7.267],label='JANAF vap') 

axes[ai,aj].plot([3130,3130],[3.554,7.460],'kD',label='C16 1-bar boiling pt.')

axes[ai,aj].set_xlabel('Temperature (K)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_title('1-bar isobar',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,10)
axes[ai,aj].set_xlim(0,4000)


#----
# density-temperature at 1 bar
ai=0
aj=1

axes[ai,aj].plot(NewEOS.onebar.T,NewEOS.onebar.rho,'+',color='blue',label='ANEOS 1 bar')

axes[ai,aj].set_ylabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_xlabel('Temperature (K)',size=sfont)
axes[ai,aj].set_title('1-bar isobar',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(1,4)
axes[ai,aj].set_xlim(0,4000)

# plot data

axes[ai,aj].plot(1999.,60.083/27.30,'s',color='red',label='Bourova Richet 1998')
axes[ai,aj].plot(1700.,60.083/23.45,'d',color='orange',label='Bourova Richet 1998')
axes[ai,aj].plot(1673.,60.083/23.51,'d',color='orange',label='')

#----
# density-pressure at 298 K
ai=1
aj=0

iT0REF = np.where(NewEOS.T == T0REF)[0]
axes[ai,aj].plot(NewEOS.rho,NewEOS.P[iT0REF[0],:],'-',color='blue',label='ANEOS 298-K isotherm')
iT0REF = np.where(NewEOS.T == 2000)[0]
axes[ai,aj].plot(NewEOS.rho,NewEOS.P[iT0REF[0],:],'-.',color='red',label='ANEOS 2000-K isotherm')
iT0REF = np.where(NewEOS.T > 2600)[0]
axes[ai,aj].plot(NewEOS.rho,NewEOS.P[iT0REF[0],:],'-.',color='cyan',label='ANEOS 2600-K isotherm')
iT0REF = np.where(NewEOS.T == T0REF)[0]

# Finkelstein et al. 2014 Forsterite isotherm
K0T = 130.0 # +-9 GPa
K0Tp = 4.12 # +- 0.07
R0F = 3.22 # g/cm3
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R0F,7./3.)-np.power(NewEOS.rho/R0F,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R0F,2./3.)-1.))
K0T = 130.0-9. # +-9 GPa
K0Tp = 4.12 # +- 0.07
BM_pressure_low = 1.5*K0T*(np.power(NewEOS.rho/R0F,7./3.)-np.power(NewEOS.rho/R0F,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R0F,2./3.)-1.))
K0T = 130.0+9. # +-9 GPa
K0Tp = 4.12 # +- 0.07
BM_pressure_high = 1.5*K0T*(np.power(NewEOS.rho/R0F,7./3.)-np.power(NewEOS.rho/R0F,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R0F,2./3.)-1.))
igood = np.where(BM_pressure < 48.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='green',label='Fo 298-K isotherm\nFinkelstein et al. 2014')
igood = np.where(BM_pressure_low < 48.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure_low[igood],'--',color='green')
igood = np.where(BM_pressure_high < 48.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure_high[igood],'--',color='green')


# Jacobsen et al. 2008 MgO isotherm
K0T = 159.6 # +-0.6 GPa room temp isotherm
K0Tp = 3.74 # +-0.03 
R00 = (24.305+15.999)*4.*1.6726/(74.697*1.E-3)/1.E3 # kg/m3->g/cm3 
print(R00)
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R00,7./3.)-np.power(NewEOS.rho/R00,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R00,2./3.)-1.))
igood = np.where((BM_pressure > 25.)&(BM_pressure < 125.))
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='brown',label='MgO 298-K isotherm\nJacobsen et al. 2008')

# Vanpeteghem et al. 2006 silicate perovskite isotherm
K0T = 253. # +-1 GPa room temp isotherm
K0Tp = 4.0 
R00 = (24.305+28.085+3.*15.999)*4.*1.6726/(162.51*1.E-3)/1.E3 # kg/m3->g/cm3 
print(R00)
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R00,7./3.)-np.power(NewEOS.rho/R00,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R00,2./3.)-1.))
igood = np.where((BM_pressure > 25.)&(BM_pressure < 125.))
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='orange',label='MgSiO$_3$ Pv 298-K isotherm\nVanpeteghem et al. 2006')

# Hosono et al. nat Geo hard sphere model isotherm 2000 K
K0T = 25. # GPa
K0Tp = 5.2
R00 = 2.65
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R00,7./3.)-np.power(NewEOS.rho/R00,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R00,2./3.)-1.))
igood = np.where(BM_pressure < 40.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='red',label='Karato 2000K hard sphere isotherm')

axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_title('298-K isotherm',size=sfont)
axes[ai,aj].set_xlim(2.5,6)
axes[ai,aj].set_ylim(0,140)

#----
if 0:
    # density-pressure at 298 K
    ai=1
    aj=1

    axes[ai,aj].plot(mantle.density,mantle.pressure,'-',color='blue',label='ANEOS '+str(np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)))+' K isentrope',markersize=10)
    axes[ai,aj].plot(PREM.density/1.E3,PREM.pressure,'-',color='orange',label='PREM')

    mantle = isentrope_class() # code in eostable.py
    Tisentrope = 1000. #K
    Smantle = np.interp(Tisentrope,NewEOS.onebar.T,NewEOS.onebar.S)

    # loop across all densities and extract the values for the requested isentrope
    for i in range(0,NewEOS.ND):
        ind = np.where((NewEOS.S[:,i] > 0))[0]
        interpfunction = interpolate.interp1d(NewEOS.S[ind,i],NewEOS.P[ind,i]) # MJ/K/kg, GPa
        mantle.pressure = np.append(mantle.pressure,interpfunction(Smantle)) # GPa
        interpfunction = interpolate.interp1d(NewEOS.S[ind,i],NewEOS.T[ind]) # MJ/K/kg, GPa
        mantle.temperature = np.append(mantle.temperature,interpfunction(Smantle)) # GPa
    mantle.density = NewEOS.rho # g/cm3

    print('Calculated isentrope with ',np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)),' K potential temperature')
    axes[ai,aj].plot(mantle.density,mantle.pressure,'-.',color='blue',label='ANEOS '+str(np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)))+' K isentrope',markersize=10)

    axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
    axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
    axes[ai,aj].tick_params(labelsize=sfont)
    axes[ai,aj].set_title('PREM and '+str(np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)))+' K ANEOS model isentrope',size=sfont)
    axes[ai,aj].set_xlim(2,6)
    axes[ai,aj].set_ylim(0,140)

# don't show a plot in lower right
#axes[1,1].axis("off")


axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)
3.6099158145574792
4.132848346563289
/var/folders/4z/xws_ywgj78sc494vpwymn9bw0000gn/T/ipykernel_87866/2785270329.py:159: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  axes[1,1].legend(fontsize=sfont-5)
Out[19]:
<matplotlib.legend.Legend at 0x163cda710>
No description has been provided for this image

ANEOS Gruneisen parameter and Theta parameter for thermal model¶

The Gruneisen parameter in forsterite liquid has been determined experimentally from shock wave data (Thomas & Asimow 2013; Davies et al. in prep.). Shockwave data derive the Grueneisen parameter in reference to a known state: \begin{equation} \gamma = \frac{1}{\rho}\frac{(P-P_{\rm ref})}{(E-E_{\rm ref})} \bigg\rvert _V \end{equation}

The ANEOS Gruneisen model is split into two parts. For $\rho \ge \rho_{0}$, \begin{equation} \gamma = \left( \frac{\gamma_0 \rho_0}{\rho} + C_{24} \left( 1-\frac{\rho_0}{\rho} \right)^2 \right) (1-C_{60}) + \left( \gamma_0 + (C_{24}-\gamma_0) \left( 1- \frac{\rho_0}{\rho} \right)^2 \right) C_{60}. \end{equation} For $\rho < \rho_{0}$, \begin{equation} \gamma = C_{16} \rho^2 + C_{17}\rho + 1 + C_{61}. \end{equation}

$C_{16}$ and $C_{17}$ are calculated so that $\gamma$ and $d \gamma / d \rho$ are continuous at $\rho=\rho_0$. The model asymptotes to the Thomas-Fermi limit ($\gamma=2/3$) when $C_{24}=2/3$.

The functional form in ANEOS does not provide a good fit to the data on forsterite. As a result, the model Gruneisen parameter does not extrapolate to the Thomas-Fermi limit at high pressures.

In [20]:
# plot the gamma and isotherms

#NewEOS.R0REF=4.0
#NewEOS.gamma0=1.1

# ANEOS Formulation for Gamma and theta
theta0=abs(NewEOS.theta0)
rhoarr = (1+np.arange(1000))/1000.*NewEOS.R0REF*10.
gammaarr = (NewEOS.gamma0*NewEOS.R0REF/rhoarr
           )+NewEOS.C24*np.power(1.0-NewEOS.R0REF/rhoarr,2.)*(1.-NewEOS.C60)+(NewEOS.gamma0 + (NewEOS.C24-NewEOS.gamma0)*np.power(1.-NewEOS.R0REF/rhoarr,2.))*NewEOS.C60
s3 = NewEOS.C61+1.
C14 = theta0 * np.exp(1.5*s3-(2.-0.5*NewEOS.C60)*NewEOS.gamma0)/(np.power(NewEOS.R0REF,s3))
C16 = (s3-(2.-NewEOS.C60)*NewEOS.gamma0)/(NewEOS.R0REF*NewEOS.R0REF)
C17 = ((3.-NewEOS.C60)*NewEOS.gamma0-2.*s3)/NewEOS.R0REF
gammaarrlow = C16*rhoarr*rhoarr+C17*rhoarr+s3
thetaarr = theta0*np.exp((1.-NewEOS.C60)*(1.-NewEOS.R0REF/rhoarr)*NewEOS.gamma0-0.5*(NewEOS.R0REF/rhoarr)*(3.-(NewEOS.R0REF/rhoarr)*(4-(NewEOS.R0REF/rhoarr))))*np.power((rhoarr/NewEOS.R0REF),NewEOS.C24)
thetalowarr = C14*np.exp(rhoarr*(C17+0.5*C16*rhoarr))*np.power(rhoarr,s3)


#-------- plots
sfont = 15
fig, axes = plt.subplots(1, 2, figsize=(14,8))
plt.subplots_adjust(wspace=0.4)


# -----
ai=0
indlow = np.where(rhoarr <= NewEOS.R0REF)
indhigh = np.where(rhoarr >= NewEOS.R0REF)
axes[ai].plot(rhoarr[indhigh[0]],gammaarr[indhigh[0]],label="ANEOS Gamma-high",color='black')
axes[ai].plot(rhoarr[indlow[0]], gammaarrlow[indlow[0]],label="ANEOS Gamma-low",color='black',ls='--')

axes[ai].plot(NewEOS.R0REF,NewEOS.gamma0,'+',label="ANEOS Gamma0",color='black')

axes[ai].set_ylim(0,2.5)
axes[ai].set_xlim(0,8)

axes[ai].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai].set_ylabel('Gamma (-)',size=sfont)
axes[ai].tick_params(labelsize=sfont)

# solid forsterite at STP
#axes[ai].plot([3.22],[1.29],'s',color='red',label='Fo solid STP')

# solid alpha-quartz at STP Rammo & Farid Powder Diffraction 1994 
axes[ai].plot([2.65],[0.641],'s',color='red',label='a-qtz solid STP')

# Thomas and Asimow 2013 liquid forsterite shock data
gamma_0=0.396
q=-2.02
q_u=1.03
Rho_0=2597. # kg/m^3
asi_rhot = np.array([3739.79122056, 3335.90237636,4691.02371542 ,4088.40354767, 4347.15217391 ,4681.1865942 ])
##asi_rhot=np.array([3340.0, 3740.0, 4100.0, 4680.0, 5290.0, 5420.0]) # BAD TABULATED VALUES in TA13
asi_gammat=gamma_0*np.power(Rho_0/asi_rhot,q)

asi_rho = np.append([2597.],asi_rhot)/1.e3
asi_gamma = np.append([0.396],asi_gammat)

#print(asi_rho)
#print(asi_gamma)
#[3739.79122056 3335.90237636 4691.02371542 4088.40354767 4347.15217391 4681.1865942 ]
#[0.82720568 0.65667847 1.30744208 0.99037675 1.12107749 1.30190973]

# Erik's error bars monte carlo'd for the asimow data
#gamma [ 0.82271901  0.65485213  1.29584059  0.97276836  1.11171104  1.29195917]
asi_gamma_sig_low=np.asarray([0., 0.25715384 , 0.14948592 , 0.5898531,   0.36154228 , 0.45607297 , 0.58790934])
asi_gamma_sig_high=np.asarray([0., 0.37855329,  0.19525818 , 1.09577441,  0.6137275 ,  0.78965792 , 1.08676142])

#axes[ai].plot(asi_rho,asi_gamma,'o',color='orange',label='Forsterite liq.\nTA13 corrected')
#for i in range(0,len(asi_gamma)):
#    axes[ai].plot([asi_rho[i],asi_rho[i]],[asi_gamma[i]-asi_gamma_sig_low[i],asi_gamma[i]+asi_gamma_sig_high[i]],linestyle='-',color='orange')

# these values from the submitted paper
#axes[ai].errorbar(SR_rho,SR_gamma,yerr=SR_gamma_e, xerr=SR_rho_e, fmt='g+', label ='Forsterite liq. Z data\nDavies et al. 2019')

gamma_0=0.5
Rho_0=2.65
Karato_gammar=np.arange(.7,1,.01)
q=-2
Karato_gammat=gamma_0*np.power(Karato_gammar,q)
#axes[ai].plot(Rho_0/(Karato_gammar),Karato_gammat,label='Karato gamma model')
#print(Rho_0/Karato_gammar)
#print(Karato_gammat)

axes[ai].plot(qtz_hug_McCoy16.rho,qtz_hug_McCoy16.gamma,'+',color='red',label=qtz_hug_McCoy16.label)
axes[ai].plot(karkigcvorig.rho,karkigcvorig.gamma,'D',label='Karki 2007')
# -----
ai=1
axes[ai].plot(rhoarr[indhigh[0]],thetaarr[indhigh[0]],label="ANEOS Theta-high",color='black')
axes[ai].plot(rhoarr[indlow[0]],thetalowarr[indlow[0]],label="ANEOS Theta-low",color='black',ls='--')
axes[ai].plot(NewEOS.R0REF,theta0,'+',label="Theta0",color='black')

axes[ai].set_ylim(0,10000)
axes[ai].set_xlim(0,8)

axes[ai].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai].set_ylabel('Theta (K)',size=sfont)
axes[ai].tick_params(labelsize=sfont)


axes[0].legend(fontsize=sfont-5)
axes[1].legend(fontsize=sfont-5)
Out[20]:
<matplotlib.legend.Legend at 0x318db5e50>
No description has been provided for this image

Specific Heat Capacity and Sounds Speeds¶

The calculated heat capacities on the forsterite Hugoniot in the liquid region, from Root et al. 2018, are similar to the 1-bar value derived from Thomas and Asimow 2013. Both are larger than the 3nR high temperature value for the solid (Gillet et al. 1991).

In [21]:
sfont = 15
fig, axes = plt.subplots(2, 2, figsize=(14,14))
plt.subplots_adjust(wspace=0.4)


#-----------------------------
# Pick isotherms and set color scale
ixtemparr = [298.,1000.,2050.,2300.,5000.,8000.,20000.]
# set color scale same for all plots
tcolmin=0.
tcolmax=20000.

# find the index closest to the selected temperatures -- no interpolation
ixindarr=np.zeros(len(ixtemparr))
ixlabelarr = []
for i in np.arange(len(ixtemparr)):
    tmpi=int(np.round(np.interp(ixtemparr[i],NewEOS.T,np.arange(NewEOS.NT))))
    #print(tmpi,antarr[tmpi])
    ixindarr[i] = tmpi
    ixlabelarr.append(str(ixtemparr[i])+' K')


#-----------------------
ai=0
aj=0
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.cv[int(ixindarr[i]),:]*1000.,c=col,linestyle='-',label=ixlabelarr[i])

#axes[ai,aj].plot(QMDH_r1,QMDH_cv*1000.,styleqmd[0],markersize=styleqmd[2],color=styleqmd[1],label='QMD Hugoniot\n(Root et al. 2018')
#axes[ai,aj].plot(2.597,1.73736,'^',color='orange',label='1 bar 2273 K liquid Cv\n(Thomas & Asimow 20133)')
#axes[ai,aj].plot(3.,180./(NewEOS.FMW/1.E3)/1.E3,'s',color='orange',label='1 bar 2000 K solid Cv\n(Gillet et al. 1991)')
axes[ai,aj].set_ylim(1.,5.)
axes[ai,aj].set_xlim(0,10)

axes[ai,aj].set_ylabel('Specific heat capacity Cv (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-----------------------
ai=0
aj=1
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.cv[int(ixindarr[i]),:]*1000.*1000./7/8.314*0.14,c=col,linestyle='-',label=ixlabelarr[i])


#axes[ai,aj].plot(QMDH_r1,QMDH_cv*1.E6/7/8.314*0.14,styleqmd[0],markersize=styleqmd[2],color=styleqmd[1],label='QMD Hugoniot\n(Root et al. 2018')
#axes[ai,aj].plot(2.597,1737./7/8.314*0.140,'^',color='orange',label='1 bar 2273 K liquid Cv\n(Thomas & Asimow 2013)')
#axes[ai,aj].plot(3.,180./(NewEOS.FMW/1.E3)/7/8.314*0.140,'s',color='orange',label='1 bar 2000 K solid Cv\n(Gillet et al. 1991)')
axes[ai,aj].set_ylim(1.,10.)
axes[ai,aj].set_xlim(0,10)

axes[ai,aj].set_ylabel('Specific heat capacity Cv (nR)',size=sfont)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-----------------------
# plot sound speed along isotherms
ai=1
aj=0
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.P[int(ixindarr[i]),:],NewEOS.cs[int(ixindarr[i]),:]/1.e5,c=col,linestyle='-',label=ixlabelarr[i])

#axes[ai,aj].plot(1.e-4,2.674,'d',color='orange',label='Thomas & Asimow 2013 2273 K')

axes[ai,aj].set_ylim(0,10)
axes[ai,aj].set_xlim(-10,40)

axes[ai,aj].set_ylabel('Bulk Sound Speed (km/s)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-----------------------
# plot sound speed along isotherms
ai=1
aj=1
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.P[int(ixindarr[i]),:],NewEOS.cs[int(ixindarr[i]),:]/1.e5,c=col,linestyle='-',label=ixlabelarr[i])

#axes[ai,aj].plot(1.e-4,2.674,'d',color='orange',label='Thomas & Asimow 2013 2273 K')

axes[ai,aj].plot(qtz_hug_McCoy16.P,qtz_hug_McCoy16.cs,'+',color='red',label=qtz_hug_McCoy16.label)


axes[ai,aj].set_ylim(0,20)
axes[ai,aj].set_xlim(0,700)

axes[ai,aj].set_ylabel('Bulk Sound Speed (km/s)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)
Out[21]:
<matplotlib.legend.Legend at 0x318dc0dd0>
No description has been provided for this image
In [ ]:
 

Compare Silica EOS Models¶

In [22]:
# Load older ANEOS models for comparison
# Load in Melosh 2013 model
# READ IN NEW ANEOS MODEL and fill the extEOStable class object
# source in eostable.py
#------------------------------------------------------------------
ORIGMODELNAME = MODELNAME
subdir='silica-Melosh2007/'
MODELNAME = 'Silica_Melosh_2007'
# Header information must all be compatible with float format
MATID = 1.0        # MATID number
DATE = 070513.     # Date as a single 6-digit number YYMMDD
VERSION = 1.0      # ANEOS Parameters Version number
FMN = 30.          # Formula weight in atomic numbers for silica
FMW = 60.083       # Formula molecular weight (g/mol)
# The following define the default initial state for material in the 201 table
R0REF   = 2.65     # g/cm3 *** R0REF is inserted into the density array
K0REF   = 3.58E11  # dynes/cm2
T0REF   = 298.     # K -- *** T0REF is inserted into the temperature array
P0REF   = 1.E6     # dynes/cm2 -- this defines the principal Hugoniot calculated below

MeloshEOS  = etab.extEOStable() # FIRST make new empty EOS object
MeloshEOS.loadextsesame(subdir+'NEW-SESAME-EXT.TXT') # LOAD THE EXTENDED 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
MeloshEOS.loadstdsesame(subdir+'NEW-SESAME-STD.TXT') # LOAD THE STANDARD 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
MeloshEOS.MODELNAME = MODELNAME # string set above in user input
MeloshEOS.MDQ = np.zeros((MeloshEOS.NT,MeloshEOS.ND)) # makes the empty MDQ array
#print(MeloshEOS.units) # these are the default units for SESAME rho-T tables
#'Units: g/cm3, K, GPa, MJ/kg, MJ/kg, MJ/K/kg, cm/s, MJ/K/kg, KPA flag. 2D arrays are (NT,ND).'

# Add the header info to the table. This could be done during the loading. 
# if made from this notebook, these values are set in the user-input above.
# ** MAKE SURE THEY MATCH ANEOS.INPUT **
MeloshEOS.MATID   = MATID
MeloshEOS.DATE    = DATE
MeloshEOS.VERSION = VERSION
MeloshEOS.FMN     = FMN
MeloshEOS.FMW     = FMW
MeloshEOS.R0REF   = R0REF
MeloshEOS.K0REF   = K0REF
MeloshEOS.T0REF   = T0REF
MeloshEOS.P0REF   = P0REF
#
# Load the information from ANEOS.INPUT and ANEOS.OUTPUT
MeloshEOS.loadaneos(aneosinfname=subdir+'ANEOS.INPUT',aneosoutfname=subdir+'ANEOS.OUTPUT',silent=True)
#
# calculate the principal Hugoniot from the table
MeloshEOS.calchugoniot(r0=MeloshEOS.R0REF,t0=MeloshEOS.T0REF)
# Calculate the fused silica Hugoniot from the table
MeloshEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.,r1=2.65)
#
# calculate the 1-bar profile; loop over temp
MeloshEOS.onebar.T = np.zeros(MeloshEOS.NT)
MeloshEOS.onebar.S = np.zeros(MeloshEOS.NT)
MeloshEOS.onebar.rho = np.zeros(MeloshEOS.NT)
it0 = np.where(MeloshEOS.T >= MeloshEOS.T0REF)[0]
id0 = np.arange(MeloshEOS.ND)#np.where(MeloshEOS.rho >= 0.8*MeloshEOS.R0REF)[0]
for iit in range(0,MeloshEOS.NT):
    MeloshEOS.onebar.T[iit] = MeloshEOS.T[iit]
    MeloshEOS.onebar.S[iit] = np.interp(1.E-4,MeloshEOS.P[iit,id0],MeloshEOS.S[iit,id0])
    MeloshEOS.onebar.rho[iit] = np.interp(1.E-4,MeloshEOS.P[iit,id0],MeloshEOS.rho[id0])
    #print(iit,MeloshEOS.onebar.rho[iit],MeloshEOS.onebar.T[iit],MeloshEOS.onebar.S[iit]*1.E3)
rho0= 2.65
r1= 2.65
367 22 2.2 298.0 0.0001 0.0 0.0 3.67465878
In [23]:
# Display the initial conditions for the Hugoniot in ANEOS.OUTPUT    
print('\nRequested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):')
print('Initial Pressure (GPa):',NewEOS.anhug.P[0])
print('Initial Temperature (K):',NewEOS.anhug.T[0])
print('Initial Sp. Entropy (MJ/K/kg):',NewEOS.anhug.S[0])
print('Initial Sp. Internal Energy (MJ/kg):',NewEOS.anhug.U[0])
print('Initial Density (g/cm3):',NewEOS.anhug.rho[0])
#

#========== PLOTTING STYLES FOR LOTS OF PLOTS =================
styleqmd = ['^','red',5,'Fo QMD calc (Root et al. 2018)']
stylezdat = ['^','black',5,'Fo Z data (Root et al. 2018)']
stylegg = ['+','green',5,'Fo Gas gun data\n(Mosenfelder et al. 2007)']
styleliq = ['^','orange',7,'Liquid Fo data corrected\n(Thomas & Asimow 2013)']
styleliqhug = ['-','orange',1,'SLVT ANEOS Liquid 2273K Hugoniot']
stylezfit = ['-','green',1,'Fo Z fit (Davies et al. 2019)']
styleanhug = ['-','red',1,'SLVT ANEOS Glass-tension Hugoniot']
stylestphug = ['-','green',3,'SLVT ANEOS Qtz Hugoniot']
styleanphase = ['-','black',1,'SLVT ANEOS Phase Boundaries']
stylephasept = ['o','black',6]

styleanFS = ['--','blue',3,'SLVT ANEOS FS Hugoniot']

styleqtzhugmelosh = ['-','goldenrod',2,'M07 ANEOS Qtz Hugoniot']
stylefshugmelosh = [':','goldenrod',2,'M07 ANEOS FS Hugoniot']
styleanphasemelosh = ['-','goldenrod',1,'M07 ANEOS Phase Boundaries']


#======================LOTS OF PLOTS============================
# plot the Hugoniots for comparison
sfont = 15
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(14,22))
plt.subplots_adjust(wspace=0.4)

#-------------
# Pressure-Density Data vs. ANEOS
ai=0
aj=0
#axes[ai,aj].plot(NewEOS.anhug.rho,NewEOS.anhug.P,styleanhug[0],linewidth=styleanhug[2],color=styleanhug[1],label=styleanhug[3])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.rho,MeloshEOS.hug.P,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.rho,MeloshEOS.hugo.P,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])
     
if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',linewidth=5,label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,':',color='green',linewidth=5,label=fs_hug_RTK19.label)

if 1:
    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].rho,shockdat_list[K13qtz['id']].P,fmt='gd',markersize=3,\
             xerr=shockdat_list[K13qtz['id']].rho_std,yerr=shockdat_list[K13qtz['id']].P_std,label=K13qtz['lab'])
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)

axes[ai,aj].set_ylim(0,1200)
axes[ai,aj].set_xlim(4.,8)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure - Temperature Data vs. ANEOS
ai=0
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T/1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T/1.e3,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.P,MeloshEOS.hug.T/1.e3,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.P,MeloshEOS.hugo.T/1.e3,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T/1.e3,color='green',linewidth=5,label='Qtz Brygoo et al. 2015')
    axes[ai,aj].plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,':',color='green',linewidth=5,label=fs_hug_Hicks06.label)

    axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T/1.E3,marker='+',color='green',label='Qtz Lyzenga & Ahrens 1980')
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T/1.E3,color='black',label=stv_mc_MM15.label)
if 1:
    axes[ai,aj].errorbar(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,fmt='b+:',\
             yerr=fs_hug_Hicks06.T_std/1.e3,label='H06 FS')
    axes[ai,aj].plot(qtz_hug_Hicks06.P,qtz_hug_Hicks06.T/1.e3,'g+:',label='H06 qtz')

axes[ai,aj].set_xlim(0,1200)
axes[ai,aj].set_ylim(0,70)
axes[ai,aj].set_ylabel('Temperature (1000 K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Particle Velocity vs. Shock Velocity; Data vs. ANEOS
ai=1
aj=0

axes[ai,aj].plot(NewEOS.hug.up,NewEOS.hug.us,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.up,NewEOS.hugo.us,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.up,MeloshEOS.hug.us,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.up,MeloshEOS.hugo.us,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.up,qtz_hug_W62.us,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.up,qtz_hug_L3.us,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.up,qtz_hug_KD.us,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.up,fs_hug_RTK19.us,':',color='green',label=fs_hug_RTK19.label)

if 1:
    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].up,shockdat_list[K13qtz['id']].us,fmt='gd',markersize=3,\
             xerr=shockdat_list[K13qtz['id']].up_std,yerr=shockdat_list[K13qtz['id']].us_std,label=K13qtz['lab'])
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.up[ind],fs_hug_RTK19.us[ind],'b--',label='R19 FS',linewidth=1)


axes[ai,aj].set_xlim(0,15)
axes[ai,aj].set_ylim(0,25)
axes[ai,aj].set_xlabel('Particle Velocity (km/s)',size=sfont)
axes[ai,aj].set_ylabel('Shock Velocity (km/s)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Sp. Entropy vs Pressure; Calculated from experiments vs. ANEOS
ai=1
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.S*1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.S*1.e3,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.P,MeloshEOS.hug.S*1.e3,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.P,MeloshEOS.hugo.S*1.e3,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    # Kraus entropy fit only good to 800 GPa
    ind = np.where(qtz_hug_KD.P <= 800.)
    axes[ai,aj].plot(qtz_hug_KD.P[ind],qtz_hug_KD.S[ind],color='green',label='Qtz Kraus et al. 2012')
    Kraus_serr = (127.3+0.328*qtz_hug_KD.P[ind]-0.000193*np.power(qtz_hug_KD.P[ind],2))/1.e3 # kJ/K/kg
    axes[ai,aj].fill_between(qtz_hug_KD.P[ind],(qtz_hug_KD.S[ind]-Kraus_serr),(qtz_hug_KD.S[ind]+Kraus_serr),color=stylezfit[1],alpha=0.2)

# silicadb
if 1:
    entropyshift = 220.
    penter = 111.
    kqtzparr = np.arange(100,500) # GPa
    kqtzsarr = 1793+22.4*kqtzparr-0.0426*np.power(kqtzparr,2)+4.319e-05*np.power(kqtzparr,3)-1.742e-8*np.power(kqtzparr,4)
    kqtzserrarr = 127.3+0.328*kqtzparr-0.000193*np.power(kqtzparr,2)
    ind = np.where(kqtzparr >= penter)[0]
    axes[ai,aj].plot(kqtzparr[ind],kqtzsarr[ind]/1.e3,'g-',linewidth=1,label='K12 qtz Hugoniot-liq region')
    axes[ai,aj].plot(kqtzparr[ind],(kqtzsarr[ind]-entropyshift)/1.e3,'g-',linewidth=3,label='K12 Qtz Hugoniot shift')
    axes[ai,aj].fill_between(kqtzparr[ind],(kqtzsarr[ind]-entropyshift-kqtzserrarr[ind])/1.e3,(kqtzsarr[ind]-entropyshift+kqtzserrarr[ind])/1.e3,color=stylezfit[1],alpha=0.2)
   
axes[ai,aj].set_ylim(0,8.)
axes[ai,aj].set_xlim(0,1000)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   LOW PRESSURE + LIQUID HUGONIOT
ai=2
aj=0
# ANEOS REFERENCE POINT (STP) Hugoniots
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.rho,MeloshEOS.hug.P,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.rho,MeloshEOS.hugo.P,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,':',color='green',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].rho,shockdat_list[L83qtzhug['id']].P,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].P_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].rho,shockdat_list[W62qtz['id']].P,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[A23qtz['id']].rho,shockdat_list[A23qtz['id']].P,fmt='gD',mfc='white',label=A23qtz['lab'],\
    #            xerr=shockdat_list[A23qtz['id']].rho_std,yerr=shockdat_list[A23qtz['id']].P_std)
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].rho,shockdat_list[T20qtz['id']].P,fmt='<',mfc='white',color='grey',\
                 xerr = shockdat_list[T20qtz['id']].rho_std, yerr=shockdat_list[T20qtz['id']].P_std, label=T20qtz['lab'])
    #print(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us)
    #res = np.polyfit(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,1)
    #print(res)
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].rho,shockdat_list[L04qtz['id']].P,fmt='gv',mfc='white',\
             xerr=shockdat_list[L04qtz['id']].rho_std,yerr=shockdat_list[L04qtz['id']].P_std,label=L04qtz['lab'])

    ind = np.where(qtz_hug_KD.P>112)[0]
    axes[ai,aj].plot(qtz_hug_KD.rho[ind],qtz_hug_KD.P[ind],'g-',linewidth=3,label='KD90 liq fit')
    axes[ai,aj].errorbar(qtzdf['rho'][indqtz],qtzdf['P'][indqtz],fmt='+',color='grey',label='M80 qtz')

axes[ai,aj].set_ylim(-5,200)
axes[ai,aj].set_xlim(2.5,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   EXTRAPOLATED HIGH PRESSURE
ai=2
aj=1
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.rho,MeloshEOS.hug.P,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.rho,MeloshEOS.hugo.P,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,':',color='green',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[W62fs['id']].rho,shockdat_list[W62fs['id']].P,fmt='bx',label=W62fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[L83fshug['id']].rho,shockdat_list[L83fshug['id']].P,fmt='bo',label=L83fshug['lab'])
    axes[ai,aj].errorbar(shockdat_list[B19fs['id']].rho,shockdat_list[B19fs['id']].P,fmt='b>',xerr=shockdat_list[B19fs['id']].rho_std,\
                 yerr=shockdat_list[B19fs['id']].P_std,mfc='white',label=B19fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[T18fs['id']].rho,shockdat_list[T18fs['id']].P,fmt='b<',mfc='white',\
                 xerr=shockdat_list[T18fs['id']].rho_std,yerr=shockdat_list[T18fs['id']].P_std,label=T18fs['lab'])

    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho)/1000.,shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'])
    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho),shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'],\
    #            xerr=shockdat_list[A23fs['id']].rho_std,yerr=shockdat_list[A23fs['id']].P_std)
    axes[ai,aj].errorbar(fsdf['rho'][indfs],fsdf['P'][indfs],fmt='+',color='grey',label='M80 FS excluded')
    indhigh = np.where(fsdf['P'][indfs] > 74)[0]
    axes[ai,aj].errorbar(fsdf['rho'][indfs][indhigh],fsdf['P'][indfs][indhigh],fmt='b+',label='M80 FS included')

    #axes[ai,aj].plot(fitfslrarr,fitfslparr,'b-',linewidth=2,label='cubic fit liquid')
    #axes[ai,aj].plot(fitfssrarr,fitfssparr,'b-',label='linear fit solid')
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)
    axes[ai,aj].errorbar(shockdat_list[R19fs['id']].rho,shockdat_list[R19fs['id']].P,fmt='bd',\
             xerr=shockdat_list[R19fs['id']].rho_std,yerr=shockdat_list[R19fs['id']].P_std,label=R19fs['lab'])

axes[ai,aj].set_ylim(0,300)
axes[ai,aj].set_xlim(3.,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)





#------------------------------
# entropy-temperature
ai=3
aj=0
axes[ai,aj].plot(NewEOS.vc.Sl*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Sv*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.S*1.E3,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Ss*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.S*1.E3,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.S[0]*1.E3,NewEOS.hug.T[0],'o',color=stylestphug[1])


axes[ai,aj].plot(MeloshEOS.vc.Sl*1.E3,MeloshEOS.vc.T,styleanphasemelosh[0],color=styleanphasemelosh[1],linewidth=styleanphasemelosh[2],label=styleanphasemelosh[3])
axes[ai,aj].plot(MeloshEOS.vc.Sv*1.E3,MeloshEOS.vc.T,styleanphasemelosh[0],color=styleanphasemelosh[1],linewidth=styleanphasemelosh[2])
axes[ai,aj].plot(MeloshEOS.cp.S*1.E3,MeloshEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])


if 1:
    #plot Connolly calcs
    b = np.argsort(C16highs)
    axes[ai,aj].plot(C16highs[b]/1.e3,C16hight[b],'-',color='black',label='C16 HSC-vWRK vapor curve')
    axes[ai,aj].plot(5.382,6074,'o',color='black',markersize=5)
    axes[ai,aj].plot(C16lows/1.e3,C16lowt,'--',color='black',label='C16 HSC-CSRK vapor curve')
    axes[ai,aj].plot(5.322,5548,'o',color='black',markersize=5)
    #plt.plot(C16sils,C16silt,'-',color='orangered',label='C16 SiO(l)-vWRK')
    #plt.plot(5403,5675,'o',color='orangered')
    # plot boiling point at 1 bar
    # note Kraus S_vap too small
    axes[ai,aj].plot([3.554,7.460],[3130,3130],'kD',label='C16 1-bar boiling pt.')

axes[ai,aj].set_ylim(0,1.5*NewEOS.cp.T)
axes[ai,aj].set_xlim(0,14)
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Temperature [K]',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-7)


#-------
# pressure - temperature melt curve
ai=3
aj=1
#axes[ai,aj].plot(NewEOS.vc.Pl,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
#axes[ai,aj].plot(NewEOS.vc.Pv,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
#axes[ai,aj].plot(NewEOS.cp.P,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.P[0],NewEOS.hug.T[0],'o',color=stylestphug[1])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T,'--',linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

if 0:
    ind = np.where(stsmcsmooth.P <=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m-',linewidth=3,label='Proposed melt curve')
    ind = np.where(stsmcsmooth.P >=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m:',linewidth=3)#,label='Proposed melt curve')
    # Millot stv melt curve crossing
    axes[ai,aj].fill_between([470.,530.],[8000.,8000.],[8600.,8600.],alpha=0.5,color='blue',label='MM15')

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,9000)
axes[ai,aj].set_xlim(-5,600)
axes[ai,aj].legend(fontsize=sfont-5)




plt.savefig(ORIGMODELNAME+'-vs-Melosh-2007.pdf')
Requested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):
Initial Pressure (GPa): 0.0001
Initial Temperature (K): 298.0
Initial Sp. Entropy (MJ/K/kg): 0.00017886
Initial Sp. Internal Energy (MJ/kg): 0.043369
Initial Density (g/cm3): 2.2
No description has been provided for this image

Check the grid spacing compared to the phase boundaries¶

Colors correspond to temperature. Black lines are ANEOS phase boundaries.

Note that ANEOS by default puts a tension region in for the solid that extends to rhomin (V17) in the input file.

If the tension region is kept, it is visible as a kink connects to the low entropy extension of the melt curve in the lower left plot of density vs. specific entropy.

The lack of points below 1000 K in the vapor dome (seen in temperature vs. specific entropy upper left) corresponds to when the sublimation pressure reaches 1.E-30. At this point, ANEOS bottoms out the pressures at this value and fixes the other state variables (it looks like to constant values; will investigate this region more closely later).

The ANEOS triple point convergence is not perfect, and there are kinks in the tabulated phase boundaries right at the triple point.

In [24]:
# check the gridding of the table around the vapor curve

# Plot aneos vapor curves
sfont = 15
fig, axes = plt.subplots(nrows=5, ncols=2, figsize=(14,26))
plt.subplots_adjust(wspace=0.4)
tcolmin=0.
tcolmax=1.2*NewEOS.cp.T
tpltmax=1.2*NewEOS.cp.T

tskip=1
#-------------
ai=0
aj=0
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Temperature [K]',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

axes[ai,aj].set_ylim(0,tpltmax)
axes[ai,aj].set_xlim(0,25)

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.e3,np.zeros(NewEOS.ND)+NewEOS.T[it],'o',markersize=2,color=col)

axes[ai,aj].plot(NewEOS.vc.Sl*1.e3,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.Sv*1.e3,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.S*1.e3,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.T,color="black")


#-------------
ai=0
aj=1
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,color=col)
axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-25,1000)
axes[ai,aj].set_xlim(1.e-20,10.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=0


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.E3,NewEOS.P[it,:],'o',markersize=2,color=col)
axes[ai,aj].semilogy(NewEOS.vc.Sl*1.e3,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.Sv*1.e3,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.S*1.e3,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-10,1000)
axes[ai,aj].set_xlim(0,15)
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=1
aj=1


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,color=col)

axes[ai,aj].semilogy(NewEOS.vc.Ul,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.Uv,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.U,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-10,1000)
axes[ai,aj].set_xlim(0,25)
axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)



#-------------
ai=2
aj=0

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,np.zeros(NewEOS.ND)+NewEOS.T[it],'o',markersize=2,color=col)
axes[ai,aj].plot(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")

axes[ai,aj].set_ylim(0.,20000)
axes[ai,aj].set_xlim(0,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=2
aj=1

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].loglog(NewEOS.rho,np.zeros(NewEOS.ND)+NewEOS.T[it],'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")

axes[ai,aj].set_ylim(10.,1.E8)
axes[ai,aj].set_xlim(1.E-10,100.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=3
aj=0

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].plot(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")


axes[ai,aj].set_ylim(0.,10)
axes[ai,aj].set_xlim(0.,15.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=3
aj=1

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].semilogx(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")

axes[ai,aj].set_ylim(0.,20)
axes[ai,aj].set_xlim(1.E-15,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)




#-------------
ai=4
aj=1

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].semilogx(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")

axes[ai,aj].set_ylim(0.,30)
axes[ai,aj].set_xlim(1.E-25,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=4
aj=0

tcolmax=40000.
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].semilogx(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")

axes[ai,aj].set_ylim(0.,200)
axes[ai,aj].set_xlim(1.E-25,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#axes[0,0].legend(fontsize=sfont-5)
#axes[0,1].legend(fontsize=sfont-5)
#axes[1,0].legend(fontsize=sfont-5)
#axes[1,1].legend(fontsize=sfont-5)
No description has been provided for this image

ANEOS KPA FLAG¶

                            TABLE          ANEOS
 KPAQQ=STATE INDICATOR      =1, 1p    =1, 1p    (eos without melt)
                            =2, 2p lv =2, 2p liquid/solid plus vapor
                                      =4, 1p solid  (eos with melt)
                                      =5, 2p melt   (eos with melt)
                                      =6, 1p liquid (eos with melt)
                            =-1 bad value of temperature
                            =-2 bad value of density
                            =-3 bad value of material number

To resolve a melt curve, there should be sufficient numbers of points within the melt curve (KPA=5, red points) so that the melt curve is continuously sampled in the grid.

The pink region is where the vapor pressure is bottomed out at 1.E-30 GPa or the pressures are negative (if there is a tension region.)

Every 10th isotherm is plotted.

In [25]:
# check the ANEOS quality flag for errors in the EOS surface
uniqkarr = np.unique(NewEOS.KPA)
print("Unique values in the KPA flag array:", uniqkarr)
print("Value: # elements with this value")
#get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x == y]

for i in uniqkarr:
    ix = get_indexes(i,np.reshape(NewEOS.KPA,NewEOS.ND*NewEOS.NT))
    print(i,': ',len(ix))


    
# Plot aneos vapor curves
sfont = 15
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14,14))
plt.subplots_adjust(wspace=0.4)
kcolmin=1
kcolmax=7

tskip=10
#-------------


#-------------
ai=0
aj=0
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=2,color=col)

axes[ai,aj].plot(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")


axes[ai,aj].set_ylim(0.,20000)
axes[ai,aj].set_xlim(0,30.)


#-------------
ai=0
aj=1

for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].loglog(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].loglog(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=2,color=col)

axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")

axes[ai,aj].set_ylim(10.,1.E8)
axes[ai,aj].set_xlim(1.E-22,100.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=0
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].semilogy(NewEOS.rho[ix],NewEOS.P[itt,ix],'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].semilogy(NewEOS.rho[ix],NewEOS.P[itt,ix],'o',markersize=2,color=col)

axes[ai,aj].semilogy(NewEOS.vc.rl,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-6,2000)
axes[ai,aj].set_xlim(0,13.)


#-------------
ai=1
aj=1
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=2,color=col)
    tmp = np.where(NewEOS.P[itt,:] < 0.)
    if (len(tmp[0]) > 0):
        axes[ai,aj].plot(NewEOS.rho[tmp[0]],np.full(len(tmp[0]),NewEOS.T[itt]),'o',markersize=2,color='purple')
    tmp = np.where(NewEOS.P[itt,:]*1.e30 <= 2.)
    if (len(tmp[0]) > 0):
        axes[ai,aj].plot(NewEOS.rho[tmp[0]],np.full(len(tmp[0]),NewEOS.T[itt]),'o',markersize=2,color='hotpink')

axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")
    
axes[ai,aj].set_ylim(10.,1.E8)
axes[ai,aj].set_xlim(1.E-22,100.)
#fig.show()
Unique values in the KPA flag array: [2. 4. 5. 6.]
Value: # elements with this value
2.0 :  104934
4.0 :  213143
5.0 :  935
6.0 :  246665
Out[25]:
(1e-22, 100.0)
No description has been provided for this image

Check the Melt Curve in Tabulated EOS¶

Colored lines correspond to isotherms. Black lines are ANEOS phase boundaries.

Tabulations of the equation of state can lead to erroneous material descriptions across the melt curve. In the following plots, the isotherms should be flat (constant pressure) across the phase boundary (black lines). Every 10th isotherm in the table is shown in the plotted regions.

In [26]:
# CHECK MELT CURVE GRIDDING
sfont = 15
fig, axes = plt.subplots(nrows=5, ncols=2, figsize=(14,26))
plt.subplots_adjust(wspace=0.5)
tcolmin=0.
tcolmax=1.5*NewEOS.cp.T
tpltmax=1.5*NewEOS.cp.T

tskip=10
#-------------
ai=0
aj=0
axes[ai,aj].plot([],[],' ',label='Low pressure melt curve')

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(.1,60)
axes[ai,aj].set_xlim(2.4,4.5)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=0
aj=1
axes[ai,aj].plot([],[],' ',label='High pressure melt curve')

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(40.,150)
axes[ai,aj].set_xlim(3.9,5.4)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=0
axes[ai,aj].plot([],[],' ',label='Higher pressure melt curve')

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")


axes[ai,aj].set_ylim(180.,310)
axes[ai,aj].set_xlim(5.5,6.4)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=1
axes[ai,aj].plot([],[],' ',label='Higher pressure melt curve')

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(180.,310)
axes[ai,aj].set_xlim(7,20)

#-------------

ai=2
aj=0
axes[ai,aj].plot([],[],' ',label='Low pressure melt curve')

axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.e3,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(.1,50)
axes[ai,aj].set_xlim(2.6,4.)


#-------------

ai=2
aj=1
axes[ai,aj].plot([],[],' ',label='High pressure melt curve')

axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.e3,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(50.,150)
axes[ai,aj].set_xlim(2.6,4.2)


#-------------
ai=3
aj=0
axes[ai,aj].plot([],[],' ',label='Low pressure melt curve')

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_xlim(1,7)
axes[ai,aj].set_ylim(.1,50)

#-------------
ai=3
aj=1
axes[ai,aj].plot([],[],' ',label='High pressure melt curve')

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")
axes[ai,aj].set_ylim(50.,150)
axes[ai,aj].set_xlim(1,15)



#-------------
ai=4
aj=0
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,color=col,linestyle='-')

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(300,700)
axes[ai,aj].set_xlim(6.5,8.)

axes[ai,aj].plot([],[],' ',label='Very high pressure melt curve.')
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=4
aj=1
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")
axes[ai,aj].set_ylim(700,1200)
axes[ai,aj].set_xlim(8.25,9.5)

axes[ai,aj].plot([],[],' ',label='Very high pressures')
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#------
axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)
axes[2,0].legend(fontsize=sfont-5)
axes[2,1].legend(fontsize=sfont-5)
axes[3,0].legend(fontsize=sfont-5)
axes[3,1].legend(fontsize=sfont-5)
axes[4,0].legend(fontsize=sfont-5)
axes[4,1].legend(fontsize=sfont-5)
Out[26]:
<matplotlib.legend.Legend at 0x336595e90>
No description has been provided for this image

References¶

Andrault, D., Bolfan-Casanova, N., Bouhifd, M. A., Boujibar, A., Garbarino, G., Manthilake, G., ... & Pesce, G. (2017). Toward a coherent model for the melting behavior of the deep Earth’s mantle. Physics of the Earth and Planetary Interiors, 265, 67-81.

Asimow, P. D. (2018). Melts under extreme conditions from shock experiments, in Magmas Under Pressure, Elsevier, p. 387-418. https://doi.org/10.1016/B978-0-12-811301-1.00015-0

Asimow P. D. and Ahrens T. J. (2010) Shock compression of liquid silicates to 125 GPa: the anorthite–diopside join. J. Geophys. Res. 115, B10209. doi:10.1029/2009JB007145.

Asimow P. D., Ghiorso M. S. (1998) Algorithmic Modifications Extending MELTS to Calculate Subsolidus Phase Relations. American Mineralogist 83, 1127-1131.

Caracas, R., Stewart, S. T. (in review) No magma ocean surface after giant impacts between rocky planets.

Centolanzi, F. J., & Chapman, D. R. (1966). Vapor pressure of tektite glass and its bearing on tektite trajectories determined from aerodynamic analysis. Journal of Geophysical Research, 71(6), 1735-1749.

Chase Jr, M. W., Curnutt, J. L., Downey Jr, J. R., McDonald, R. A., Syverud, A. N., & Valenzuela, E. A. (1982). JANAF thermochemical tables, 1982 supplement. Journal of Physical and Chemical Reference Data, 11(3), 695-940. https://janaf.nist.gov

Costa, G. C., Jacobson, N. S., & Fegley Jr, B. (2017). Vaporization and thermodynamics of forsterite-rich olivine and some implications for silicate atmospheres of hot rocky exoplanets. Icarus, 289, 42-55.

de Koker, N. P., Stixrude, L., & Karki, B. B. (2008). Thermodynamics, structure, dynamics, and freezing of Mg2SiO4 liquid at high pressure. Geochimica et Cosmochimica Acta, 72(5), 1427-1441.

Davies, E.J., Carter, P.J., Root, S., Kraus, R.G., Spaulding, D.K., Stewart, S.T., Jacobsen, S.B. (2019) Silicate Melting and Vaporization during Rocky Planet Formation, J. Geophys. Res. Planets, submitted.

Duffy, T., Madhusudhan, N., & Lee, K. K. M. (2015). 2.07 Mineralogy of Super-Earth Planets. Treatise on Geophysics, Second Edition, Elsevier, Oxford.Duffy, T., Madhusudhan, N., & Lee, K. K. M. (2015). 2.07 Mineralogy of Super-Earth Planets. Treatise on Geophysics, Second Edition, Elsevier, Oxford.

Fiquet, G., Auzende, A. L., Siebert, J., Corgne, A., Bureau, H., Ozawa, H., & Garbarino, G. (2010). Melting of peridotite to 140 gigapascals. Science, 329(5998), 1516-1518.

Gillet, P., Richet, P., Guyot, F., & Fiquet, G. (1991). High‐temperature thermodynamic properties of forsterite. Journal of Geophysical Research: Solid Earth, 96(B7), 11805-11816.

Ghiorso, Mark S., and Sack, Richard O. (1995) Chemical Mass Transfer in Magmatic Processes. IV. A Revised and Internally Consistent Thermodynamic Model for the Interpolation and Extrapolation of Liquid-Solid Equilibria in Magmatic Systems at Elevated Temperatures and Pressures. Contributions to Mineralogy and Petrology, 119, 197-212.

Hosono, N., Karato, S. I., Makino, J., & Saitoh, T. R. (2019). Terrestrial magma ocean origin of the Moon. Nature Geoscience, 1, 418-423.

Jacobsen, S.D., Holl, C.M., Adams, K.A., Fischer, R.A., Martin, E.S., Bina, C.R., Lin, J.F., Prakapenka, V.B., Kubo, A. and Dera, P., 2008. Compression of single-crystal magnesium oxide to 118 GPa and a ruby pressure gauge for helium pressure media. American Mineralogist, 93(11-12), pp.1823-1828.

Jing, Z, and Karato, S. (2011) A new approach to the equation of state of silicate melts: An application of the theory of hard sphere mixtures, Geochimica et Cosmochimica Acta 75, 6780–6802.

Kraus, R.G., Stewart, S.T., Swift, D.C., Bolme, C.A., Smith, R.F., Hamel, S., Hammel, B.D., Spaulding, D.K., Hicks, D.G., Eggert, J.H. and Collins, G.W. (2012). Shock vaporization of silica and the thermodynamics of planetary impact events. Journal of Geophysical Research: Planets, 117(E9), E09009, doi:10.1029/2012JE004082.

Kumazawa, M., & Anderson, O. L. (1969). Elastic moduli, pressure derivatives, and temperature derivatives of single‐crystal olivine and single‐crystal forsterite. Journal of Geophysical Research, 74(25), 5961-5972.

Luo, S. N., Akins, J. A., Ahrens, T. J., & Asimow, P. D. (2004). Shock‐compressed MgSiO$_3$ glass, enstatite, olivine, and quartz: Optical emission, temperatures, and melting. Journal of Geophysical Research: Solid Earth, 109(B5).

Lyon, S. P., & Johnson, J. D. (1992). SESAME: The LANL equation of state database. Los Alamos National Laboratories Report LAUR-92-3407, Los Alamos, NM.

McDonough, W. and Sun, S.-s. (1995). The composition of the Earth. Chemical Geology 120, 223-253.

Marinova, M. M., Aharonson, O., & Asphaug, E. (2011). Geophysical consequences of planetary-scale impacts into a Mars-like planet. Icarus, 211(2), 960-985.

Mosenfelder, J. L., Asimow, P. D., & Ahrens, T. J. (2007). Thermodynamic properties of Mg2SiO4 liquid at ultra‐high pressures from shock measurements to 200 GPa on forsterite and wadsleyite. Journal of Geophysical Research: Solid Earth, 112(B6).

Mysen, B. O., & Kushiro, I. (1988). Condensation, evaporation, melting, and crystallization in the primitive solar nebula; experimental data in the system MgO-SiO$_2$-H$_2$ to $1.0\times10^{-9}$ bar and 1870 degrees C with variable oxygen fugacity. American Mineralogist, 73(1-2), 1-19.

Nagahara, H., Kushiro, I., & Mysen, B. O. (1994). Evaporation of olivine: Low pressure phase relations of the olivine system and its implication for the origin of chondritic components in the solar nebula. Geochimica et cosmochimica acta, 58(8), 1951-1963.

Richet, P., Leclerc, F., & Benoist, L. (1993). Melting of forsterite and spinel, with implications for the glass transition of Mg$_2$SiO$_4$ liquid. Geophysical Research Letters, 20(16), 1675-1678.

Robie, R. A., Hemingway, B. S., & Takei, H. (1982). Heat capacities and entropies of Mg$_2$SiO$_4$, Mn$_2$SiO$_4$, and Co$_2$SiO$_4$ between 5 and 380 K. American Mineralogist, 67(5-6), 470-482.

Root, S., Townsend, J.P., Davies, E., Lemke, R.W., Bliss, D.E., Fratanduono, D.E., Kraus, R.G., Millot, M., Spaulding, D.K., Shulenburger, L. and Stewart, S.T. (2018). The principal Hugoniot of forsterite to 950 GPa. Geophysical Research Letters, 45(9), 3865-3872.

Solomatova, N. V., & Caracas, R. (2019). Pressure‐induced coordination changes in a pyrolitic silicate melt from ab initio molecular dynamics simulations. Journal of Geophysical Research: Solid Earth, 124, 11232–11250. https://doi.org/10.1029/2019JB018238

Suzuki, A., & Ohanti, E. (2003). Density of peridotite melts at high pressure. Phys. Chem. Minerals 30, 449-456, doi:10.1007/s00269-003-0322-6.

Thomas, C. W., & Asimow, P. D. (2013). Direct shock compression experiments on premolten forsterite and progress toward a consistent high‐pressure equation of state for CaO‐MgO‐Al$_2$O$_3$‐SiO$_2$‐FeO liquids. Journal of Geophysical Research: Solid Earth, 118(11), 5738-5752.

Tillotson, J. H. (1962). Metallic equations of state for hypervelocity impact (No. GA-3216). General Atomics Division, General Dynamics, San Diego, CA.

Townsend, J. P., Shohet, G., & Cochrane, K. R. (2020). Liquid-vapor coexistence and critical point of Mg2SiO4 from ab initio simulations. Geophysical Research Letters, 47, e2020GL089599. https:// doi.org/10.1029/2020GL089599

Vanpeteghem, C. B., Zhao, J., Angel, R. J., Ross, N. L., & Bolfan‐Casanova, N. (2006). Crystal structure and equation of state of MgSiO3 perovskite. Geophysical Research Letters, 33(3), 2005GL024955.

Xiao, B., & Stixrude, L. (2018). Critical vaporization of MgSiO$_3$. Proceedings of the National Academy of Sciences, 115(21), 5371-5376.

Zeman, M., Holec, M., & Váchal, P. (2019). HerEOS: A framework for consistent treatment of the Equation of State in ALE hydrodynamics. Computers & Mathematics with Applications, 78(2), 483-503.

ANEOS references¶

Collins, Gareth S., and H. Jay Melosh (2014). Improvements to ANEOS for multiple phase transitions. 45th Lunar Planet. Sci. Conf. Abs. 2664.

Melosh, H. J. (2007). A hydrocode equation of state for SiO$_2$. Meteoritics & Planetary Science, 42(12), 2079-2098.

Thompson, S. L. (1990). ANEOS analytic equations of state for shock physics codes input manual. SANDIA REPORT SAND, 89-2951.

Thompson, S. L., & Lauson, H. S. (1974). Improvements in the Chart D radiation-hydrodynamic CODE III: Revised analytic equations of state (No. SC-RR--71-0714). Sandia Labs.

Thompson, S. L., Lauson, H. S., Melosh, H. J., Collins, G. S., & Stewart, S. T. (2019, November 1). M-ANEOS (Version 1.0). Zenodo. http://doi.org/10.5281/zenodo.3525030

Stewart, S. T. (2019). ANEOS Code Modification: Thermal model adjustment parameter. https://github.com/ststewart/aneos-forsterite-2019/EOS-docs/

ANEOS Input Parameters¶

EOS -1: low temperature solid model
THUG, RHUG: initial state for the Hugoniot calculated in ANEOS.OUTPUT.

V01 nelem = 1: This model has 1 elements: Si, O.
V02 model type=4: solid-liquid-gas model with ionization
V03 rho0=2.2 g/cm3: reference solid density at the reference state, in this case STP.
V04 T0=298 K: reference state is STP.
V05 P0=1.e6 dynes/cm2: reference state STP.
V06 B0=0.40E12 dynes/cm2: bulk modulus is close to alpha-quartz true value.
V07 gamma0=0.57: Grueneisen parameter.
V08 Tdebye=-1600.0 K: fitted to obtain an entropy of vaporization for silica. Negative means use full Debye model.

V09 TG model=1: based on free-volume theory, generate a derivative of the bulk modulus near 4.
V10 3*C24=6.0: Compromise value for the gamma function to span solid to liquid.
V11 Esep=1.2E11 erg/g: Zero temperature separation energy. Fitted to match the triple point pressure and critical point temperature.
V12 Tmelt=1996 K: Melting temperature at reference pressure. Here, using the value for alpha-quartz from Bourova and Richet 1998..
V13 C53=0 erg/g: Critical point adjustment parameter. Not used.
V14 C54=0. [dimless]: Critical point adjustment parameter. Not used.
V15 H0=0: Thermal conduction parameter. Not used.
V16 C41=0: Thermal conduction parameter. Not used.

V17 rhomin=0 g/cm3: minimum density for the solid. Defines the tension region in the model. The default value is 0.8*rho0 when rhomin=0.
V18=2.92 g/cm3: Solid-solid phase transition parameter. Density at onset of transition.
V19=4.30 g/cm3: Solid-solid phase transition parameter. Density at end of transition.
V20=0.4E11 dynes/cm2: Solid-solid phase transition parameter. Pctr, pressure at the center of the transition.
V21=2.85E12: Solid-solid phase transition parameter. First derivative of Pctr with respect to density. Bulk modulus of stishovite. https://doi.org/10.1016/0167-6636(93)90083-4

V22=1.E1: Solid-solid phase transition parameter. Second derivative of Pctr with respect to density.
V23 Hfusion=0.66E10 erg/g: Latent heat of fusion at the reference pressure. Fitted to provide a melt curve similar to Millot et al. 2015.
V24 rhol/rhos=0.93: Volume change on melting. Compromise value.

V25 upper=0.0: Upper limit to cold compression curve extension. Default=1. Setting to zero means use default value.
V26 lower=0.: Default=0.
V27 alpha=0.3: Liquid model parameter ($0<\alpha<1$, default=0.3).
V28 beta=0.1: Liquid model parameter ($0<\beta<1$, default=0.1). beta cannot be equal to gamma.
V29 gamma=0.3: Liquid model parameter ($0<\gamma<1$, default=0.2). beta cannot be equal to gamma.
V30 C60=0: Gamma model adjustment parameter. Default=0.
V31 C61=-0.851 Gamma model adjustment parameter to provide best fit to the critical point. Default=0.
V32 C62=0.5: Critical point adjustment parameter. Default=0. ($0

V33 Ionization model=0=Saha model.
V34-V35=0: Reactive chemistry model not used.
Melosh molecular clusters model for the critical point (Melosh fitted values for V36-V43):
V36 Natom=2: number of atoms in molecular clusters
V37 Ebind=5.0 eV: Binding energy
V38 RotDOF=2.0: Rotational degrees of freedom, 2 for diatomic molecule
V39 Rbond=1.5E-8 cm: Length of molecular bond (cm)
V40 VibDOF=1: Numer of vibrational degrees of freedom, 1 for diatomic molecule

V41 Tdebye=2000 K: Vibrational Debye temperature
V42 Mieflag=1: Flag for Mie potential (1) or Morse potential.
V43 a_exp=1.7: Power in Mie potential (1 to 2). Fitted to improve critical point
V44 $f_{cv}=1.40$: Adjust heat capacity in high temperature limit. $C_v=3f_{cv}NkT$.
V45 QCC1: 1.E-30. low density value to transition to ideal gas.
V45 QCC6: 1.E5. high temperature psi value to transition to ideal gas.
V46 to V48=0: Variables not used.

Atomic number and atom fraction:
8 0.66667
14 0.33333

Development Notes¶

This preliminary fused silica model was developed using the data from K. Amodeo's PhD thesis (manuscript in prep.).

End of File

¶